Disclaimer: The purpose of the Open Case Studies project is to demonstrate the use of various data science methods, tools, and software in the context of messy, real-world data. A given case study does not cover all aspects of the research process, is not claiming to be the most appropriate way to analyze a given dataset, and should not be used in the context of making policy decisions without external consultation from scientific experts.
Body mass index (BMI) is often used as a proxy for adiposity (the condition of having excess body fat) and is measured as an individual’s weight in kilograms (kg) or pounds (lbs) divided by the individual’s height in meters (\(m^2\)) squared.
There are different categories of weight defined using different BMI thresholds as defined by the World Health Organization (WHO), the Centers for Disease Control and Prevention, and the National Institutes of Health. For example, these are the categories defined by the WHO:
The following chart can help you identify your estimate of BMI: 
Previous work has shown higher BMI (>30) is associated with higher rates of all-causes of mortality, as well as increased rates of type 2 diabetes, cancer, heart disease, and stroke.
Recently, an article published in Nature evaluated and compared the Body-Mass Index (BMI) of populations in rural and urban communities around the world:
NCD Risk Factor Collaboration (NCD-RisC). Rising rural body-mass index is the main driver of the global obesity epidemic in adults. Nature 569, 260–264 (2019).
The article challenged the widely-held view that increased urbanization was one of the major reasons for increased global obesity rates. This view came about because many countries around the world have shown increased urbanization levels in parallel with increased obesity rates. In this article, however, the NCD-RisC argued that this might not be the case and that in fact for most regions around the world, BMI measurements are increasing in rural populations just as much, if not more so, than urban populations. Furthermore, this study suggested that obesity has particularly increased in female populations in rural regions:
“We noted a persistently higher rural BMI, especially for women.”
In this case study, we will evaluate the data reported in this article to explore regional and gender specific differences in the obesity rates around the world in 1985 and 2017. Most importantly we will test if there is a difference in obesity rates between rural and urban communities and if there has been a change in obesity over time for these regions, particularly for women. To do this we will test if there is a difference between the average BMI for each group.
While gender and sex are not actually binary, the data presented that is used in this analysis only contain data for groups of individuals described as men or women.
Our main questions are:
In this case study, we will walk you through importing data from a pdf, cleaning, wrangling and visualizing the data, and comparing two groups .
We will use well-established and commonly used packages, including stringr, tidyr, dplyr, purrr, and ggplot2 from the tidyverse. Specifically, the tidyverse is
“an opinionated collection of R packages designed for data science. All packages share an underlying philosophy and common APIs”.
Another way of putting it is that it’s a set of packages that are useful specifically for data manipulation, exploration and visualization with a common philosophy. While some students may be familiar with previous packages from the R programming language, these packages make data science in R especially efficient.
We will begin by loading the packages that we will need with a brief explanation:
library(here)
library(pdftools)
library(stringr)
library(readr)
library(dplyr)
library(tibble)
library(magrittr)
library(glue)
library(purrr)
library(tidyr)
library(ggplot2)
library(ggrepel)
library(cowplot)
library(patchwork)| Package | Use |
|---|---|
| here | to easily load and save data with relative paths |
| pdftools | to read a text from pdf into R |
| stringr | to manipulate the text data |
| readr | to manipulate the text data within the pdf into individual lines |
| dplyr | to arrange/filter/select subsets of the data |
| tibble | to create data objects that we can manipulate with dplyr/stringr/tidyr/purrr |
| magrittr | to use the %<>% piping operator |
| glue | to paste or combine character strings and data together |
| purrr | to perform functions on all columns of a tibble |
| tidyr | to convert data from ‘wide’ to ‘long’ format |
| ggplot2 | to make visualizations with multiple layers |
| ggrepel | to allow labels in figures not to overlap |
| cowplot and patchwork | to allow plots to be combined |
The first time we use a function, we will use the :: to indicate which package we are using (e.g. stringr::str_detect()). Unless we have overlapping function names, this is not necessary, but we will include it here to be informative about where the functions we will use come from.
The measurement of BMI has some limitations that are well recognized, as it does not account for the composition of body mass, the location of body fat, or the contribution of body frame size. However, BMI has been a useful health indicator for risk for many diseases and conditions particularly when combined with other risk factor information.
We will be using data located within a table of the supplementary material for the NCD-RisC paper referenced above.
This is a pdf that can be found freely available online. Here is a screenshot of the first few rows of this table:
You can see that the data contain average BMI estimates for both men and women in countries at the national level, as well as the average BMI estimates for the rural and urban areas of these countries in both 1985 and 2017.
The data within the parentheses are the 95% credible intervals (CIs) for the average BMI estimates. The authors provide these CIs as a guide to understand how likely the estimate is for the true population mean BMI. A wider range suggests that the estimate is less accurate, as there are more possible values for the true mean with credible evidence. Notice that the mean BMI values are “age-standardised,” this means that the mean values were adjusted for the different age distributions of the various countries so that the countries can be more fairly compared. For example, if one country has a population that is considerably younger, the mean BMI might be quite low (as younger individuals tend to have lower BMI values due to faster metabolic rates). One might falsely conclude that people in that country generally have lower BMI values than people in most other countries, however the lower overall BMI mean might simply be due to the fact that the country has a younger population than those other countries.
First let’s download the data from the PDF file. We use file.exists() from base R to check if the file we want to download already exists. If not, then we download it using the function utils::download.file(). We also use the here() function of the here package to look specifically in the docs directory` of the directory where our open RStudio project file is located. This helps us, as we do not need to write the full path to the document. This makes our code more reproducible for someone else, who may have a different overall file structure on their computer but downloads our case study files from GitHub and therefore has the same file structure specifically for the files included in our GitHub repository. Thus our commands about file locations still work on their local machine.
if(!file.exists(here("docs", "paper_supplement.pdf"))){
url <- "https://static-content.springer.com/esm/art%3A10.1038%2Fs41586-019-1171-x/MediaObjects/41586_2019_1171_MOESM1_ESM.pdf"
utils::download.file(url, here("docs", "paper_supplement.pdf"))
}Now that we have downloaded the PDF file, we will read it in to R using the pdftools package:
Let’s take a look at the data – the summary() function from base R helps us to look at the structure of R objects.
Length Class Mode
63 character character We can see that we have 63 elements that are character strings. You may also notice that the original PDF has 63 pages. Let’s take a look at some of these elements and compare them to the pages in the original document.
[1] "Letter https://doi.org/10.1038/s41586-019-1171-x\nRising rural body-mass index is the main driver of\nthe global obesity epidemic in adults\nNCD Risk Factor Collaboration (NCD-RisC)*\n*A list of authors and their affiliations appears in the online version of the paper.\n2 6 0 | N A T U RE | V O L 5 6 9 | 9 M A Y 2 0 1 9\n"[1] "Supplementary Information. Statistical model for estimating BMI trends by rural and urban\nplace of residence.\n 1\n" chr " 2 2\n Age-standardised mean BMI in 1985 (kg/m ) Age-standardised mean BMI in 2017 (kg/m )\n Country Sex\n National Rural Urban National Rural Urban\n Men 20.2 (17.8-22.7) 19.7 (17.2-22.2) 22.4 (20.0-25.0) 22.8 (20.3-25.3) 22.5 (20.0-25.0) 23.6 (21.0-26.1)\nAfghanistan Women 20.6 (18.4-22.8) 20.1 (17.8-22.4) 23.2 (20.9-25.4) 24.4 (23.3-25.4) 23.6 (22.5-24.8) 26.3 (25.1-27.4)\n Men 25.2 (23.9-26.5) 25.0 (23.7-26.4) 25.4 (24.0-26.7) 27.0 (26.0-27.9) 26.9 (25.9-27.9) 27.0 (26.0-28.0)\nAlbania Women 26.0 (24.1-27.9) 26.1 (24.1-28.1) 25.9 (23.9-27.8) 26.0 (24.8-27.2) 26.2 (24.8-27.5) 25.9 (24.6-27.2)\n Men 22.1 (20.8-23.3) 21.8 (20.5-23.1) 22.3 (21.0-23.6) 25.1 (24.5-25.7) 24.8 (24.1-25.4) 25.2 (24.6-25.9)\nAlgeria Women 24.0 (22.2-25.7) 23.3 (21.4-25.1) 24.8 (22.9-26.6) 27.4 (26.7-28.0) 27.0 (26.3-27.8) 27.5 (26.7-28.2)\n Men 33.7 (32.7-34.7) 32.6 (31.7-33.5) 34.0 (32.9-35.1) 34.3 (33.0-35.6) 34.6 (33.1-35.9) 34.2 (32.9-35.6)\nAmerican Samoa Women 34.3 (33.1-35.6) 34.1 (33.0-35.2) 34.4 (33.0-35.8) "| __truncated__Yes, this looks like data in a the table we want. We are using the str() function to print just a portion of the data. The str() function compactly displays the structure of arbitrary R objects. The nchar.max argument let’s us define the maximal number of characters to show. Note print(pdf_obesity[52]) would show the entire contents of the this page of the table.
Here is the PDF version again:
We can see that the output looks pretty similar to the pages of the pdf, but the spacing is a bit awkward. Note that the way the data is displayed is partially influenced by the width setting of the RStudio window.
Now, we are interested in a supplementary Table 3, which has multiple pages and includes the same header on each page. We can use that to determine what elements of our pdf_obesity character strings include our table. We will use the str_detect() function from the stringr package to search for the elements that contain text that is consistently in the header. The output of of this function will show which elements of the object (in this case pages of the pdf) include this pattern indicated as a TRUE or FALSE.
[1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[37] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[49] FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[61] TRUE FALSE FALSEWe note, the “Age-standardised mean BMI” is part of the header in the table on every page. The code above shows the pages that contain our table of interest (as TRUE values).
Now we will extract just the data for the table now and call it rural_urban. To do this we will use the str_subset() function of the stringr package.
rural_urban <- stringr::str_subset(pattern = "Age-standardised mean BMI",
string = pdf_obesity)
summary(rural_urban) Length Class Mode
10 character character Here, we see there are 10 pages worth of elements in our rural_urban object.
Let’s check the first and last page:
chr " 2 2\n Age-standardised mean BMI in 1985 (kg/m ) Age-standardised mean BMI in 2017 (kg/m )\n Country Sex\n National Rural Urban National Rural Urban\n Men 20.2 (17.8-22.7) 19.7 (17.2-22.2) 22.4 (20.0-25.0) 22.8 (20.3-25.3) 22.5 (20.0-25.0) 23.6 (21.0-26.1)\nAfghanistan Women 20.6 (18.4-22.8) 20.1 (17.8-22.4) 23.2 (20.9-25.4) 24.4 (23.3-25.4) 23.6 (22.5-24.8) 26.3 (25.1-27.4)\n Men 25.2 (23.9-26.5) 25.0 (23.7-26.4) 25.4 (24.0-26.7) 27.0 (26.0-27.9) 26.9 (25.9-27.9) 27.0 (26.0-28.0)\nAlbania Women 26.0 (24.1-27.9) 26.1 (24.1-28.1) 25.9 (23.9-27.8) 26.0 (24.8-27.2) 26.2 (24.8-27.5) 25.9 (24.6-27.2)\n Men 22.1 (20.8-23.3) 21.8 (20.5-23.1) 22.3 (21.0-23.6) 25.1 (24.5-25.7) 24.8 (24.1-25.4) 25.2 (24.6-25.9)\nAlgeria Women 24.0 (22.2-25.7) 23.3 (21.4-25.1) 24.8 (22.9-26.6) 27.4 (26.7-28.0) 27.0 (26.3-27.8) 27.5 (26.7-28.2)\n Men 33.7 (32.7-34.7) 32.6 (31.7-33.5) 34.0 (32.9-35.1) 34.3 (33.0-35.6) 34.6 (33.1-35.9) 34.2 (32.9-35.6)\nAmerican Samoa Women 34.3 (33.1-35.6) 34.1 (33.0-35.2) 34.4 (33.0-35.8) "| __truncated__Using [1] shows just the first page of the 10 pages that contain the table. Also, notice the new line "\n" values that indicate new lines.
This looks the same as the beginning… how about the end?
chr " 2 2\n Age-standardised mean BMI in 1985 (kg/m ) Age-standardised mean BMI in 2017 (kg/m )\n Country Sex\n National Rural Urban National Rural Urban\n Men 20.0 (18.6-21.4) 19.6 (18.1-21.1) 20.7 (19.2-22.0) 22.2 (21.6-22.9) 21.7 (21.0-22.4) 23.0 (22.3-23.7)\nZambia Women 21.4 (20.3-22.6) 20.7 (19.5-21.9) 22.5 (21.3-23.7) 23.9 (23.3-24.5) 22.7 (22.1-23.3) 25.5 (24.9-26.2)\n Men 21.2 (20.2-22.3) 20.9 (19.8-22.1) 22.1 (21.1-23.1) 22.3 (21.7-22.9) 21.9 (21.2-22.5) 23.3 (22.6-23.9)\nZimbabwe Women 24.6 (23.5-25.9) 24.0 (22.8-25.4) 26.5 (25.3-27.8) 25.3 (24.6-26.1) 24.5 (23.7-25.3) 27.1 (26.2-27.9)\n*The entire population live in areas classified as urban (Bermuda, Hong Kong (from 1992 onwards), Nauru and Singapore) or rural (Tokelau).\nNumbers in parentheses show 95% credible intervals.\n 60\n"
Great! Our rural_urban object looks like it contains the entire Supplementary 3 table, as both the beginning and the end include the data we expected.
At this point we have large strings now for each page of the table, but this is not very convenient to work with. Now, we will wrangle the data into a more usable form. Ideally, we would like our data to be in some sort of tabular form.
First, it would be useful to separate the data that is currently in 10 giant string chunks into individual lines or rows of the table. To do this we can use the read_lines() function of the readr package which will separate the lines based on the "\n" new line characters.
Length Class Mode
10 character character Length Class Mode
461 character character We can see that each string in the rural_urban object is now a single line from the table.
For example, line 6 shows the data for women in Afghanistan.
[1] "Afghanistan Women 20.6 (18.4-22.8) 20.1 (17.8-22.4) 23.2 (20.9-25.4) 24.4 (23.3-25.4) 23.6 (22.5-24.8) 26.3 (25.1-27.4)"We also have a lot of white-space. Let’s get rid of the excess white spaces using str_squish() also from the stringr package.
[1] "2 2"
[2] "Age-standardised mean BMI in 1985 (kg/m ) Age-standardised mean BMI in 2017 (kg/m )"
[3] "Country Sex"
[4] "National Rural Urban National Rural Urban"
[5] "Men 20.2 (17.8-22.7) 19.7 (17.2-22.2) 22.4 (20.0-25.0) 22.8 (20.3-25.3) 22.5 (20.0-25.0) 23.6 (21.0-26.1)"
[6] "Afghanistan Women 20.6 (18.4-22.8) 20.1 (17.8-22.4) 23.2 (20.9-25.4) 24.4 (23.3-25.4) 23.6 (22.5-24.8) 26.3 (25.1-27.4)"Now it is much easier to see the data.
If we look at the part of our data that contains the end of the first page of the table and the start of the second we can see that the header information is repeated, as well as a line with the page number and an empty line, and a line that says “2 2”.
[1] "Benin Women 20.6 (19.1-22.0) 20.1 (18.5-21.6) 21.8 (20.3-23.3) 24.3 (23.5-25.0) 23.2 (22.5-24.0) 25.5 (24.7-26.3)"
[2] "Men 24.3 (21.8-26.8) na* 24.3 (21.8-26.8) 26.7 (24.3-29.2) na* 26.7 (24.3-29.2)"
[3] "Bermuda Women 25.4 (22.2-28.6) na* 25.4 (22.2-28.6) 28.5 (25.3-31.6) na* 28.5 (25.3-31.6)"
[4] "Men 20.6 (19.1-22.1) 20.3 (18.7-21.8) 23.1 (21.5-24.6) 23.5 (22.8-24.3) 23.1 (22.2-23.9) 24.2 (23.3-25.2)"
[5] "Bhutan Women 20.7 (18.4-22.9) 20.3 (18.0-22.6) 23.0 (20.8-25.3) 24.6 (23.7-25.5) 23.8 (22.7-24.9) 25.9 (24.7-27.0)"
[6] "51"
[7] ""
[8] "2 2"
[9] "Age-standardised mean BMI in 1985 (kg/m ) Age-standardised mean BMI in 2017 (kg/m )"
[10] "Country Sex"
[11] "National Rural Urban National Rural Urban"
[12] "Men 23.3 (21.0-25.5) 22.9 (20.6-25.2) 23.6 (21.3-25.9) 26.1 (23.9-28.4) 25.3 (23.1-27.5) 26.5 (24.2-28.8)"
[13] "Bolivia Women 23.8 (22.3-25.3) 23.0 (21.4-24.5) 24.6 (23.1-26.1) 27.9 (26.5-29.3) 27.0 (25.6-28.5) 28.3 (26.8-29.7)"Although the header was necessary on all of the pages of the pdf version of the table, we only need that information once in our data.
So, let’s remove all the header information and the page number lines from the rural_urban object, then we will make a single line header for the beginning. One way to do this is to find all lines that include either “Women” or “Men” and only keep this data.
First, let’s see how many lines include “Women” or “Men”. We can do this by first identifying the lines that contain these patterns and then check the length of the resulting vector for the line numbers that contain these patterns. Note, the "|" indicates that we are looking for either “Women” or “Men”.
[1] 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
[19] 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
[37] 41 42 43 44 45 46 47 48 55 56 57 58 59 60 61 62 63 64
[55] 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82
[73] 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 105 106
[91] 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124
[109] 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142
[127] 143 144 145 146 147 148 155 156 157 158 159 160 161 162 163 164 165 166
[145] 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184
[163] 185 186 187 188 189 190 191 192 193 194 195 196 197 198 205 206 207 208
[181] 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226
[199] 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244
[217] 245 246 247 248 255 256 257 258 259 260 261 262 263 264 265 266 267 268
[235] 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286
[253] 287 288 289 290 291 292 293 294 295 296 297 298 305 306 307 308 309 310
[271] 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328
[289] 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346
[307] 347 348 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370
[325] 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388
[343] 389 390 391 392 393 394 395 396 397 398 405 406 407 408 409 410 411 412
[361] 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430
[379] 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448
[397] 455 456 457 458[1] 400Note that str_which() is case sensitive, so it would not work to use “women” as the pattern, and using “men” would return the lines that contain “Women” or “Yemen” etc. Try this out with pattern = "women" and pattern = "men" to see the result!
OK, so this looks correct. This includes most lines but there are gaps where the header is located. It looks like out of the original 461 lines, there are 400 lines in our table that aren’t headers.
We can check our data now using either head() from base R or glimpse() from the tibble and dplyr package. (Yup, the same function is available from both packages!) The head() function shows us the first rows or lines of the data, while the glimpse() function provides us information about the total size of the object and shows us the first line or row.
chr [1:400] "Men 20.2 (17.8-22.7) 19.7 (17.2-22.2) 22.4 (20.0-25.0) 22.8 (20.3-25.3) 22.5 (20.0-25.0) 23.6 (21.0-26.1)" ...Now our rural_urban object 400 lines instead of 461, like it did before when it contained the redundant header information.
[1] "Men 20.2 (17.8-22.7) 19.7 (17.2-22.2) 22.4 (20.0-25.0) 22.8 (20.3-25.3) 22.5 (20.0-25.0) 23.6 (21.0-26.1)"
[2] "Afghanistan Women 20.6 (18.4-22.8) 20.1 (17.8-22.4) 23.2 (20.9-25.4) 24.4 (23.3-25.4) 23.6 (22.5-24.8) 26.3 (25.1-27.4)"
[3] "Men 25.2 (23.9-26.5) 25.0 (23.7-26.4) 25.4 (24.0-26.7) 27.0 (26.0-27.9) 26.9 (25.9-27.9) 27.0 (26.0-28.0)"
[4] "Albania Women 26.0 (24.1-27.9) 26.1 (24.1-28.1) 25.9 (23.9-27.8) 26.0 (24.8-27.2) 26.2 (24.8-27.5) 25.9 (24.6-27.2)"
[5] "Men 22.1 (20.8-23.3) 21.8 (20.5-23.1) 22.3 (21.0-23.6) 25.1 (24.5-25.7) 24.8 (24.1-25.4) 25.2 (24.6-25.9)"
[6] "Algeria Women 24.0 (22.2-25.7) 23.3 (21.4-25.1) 24.8 (22.9-26.6) 27.4 (26.7-28.0) 27.0 (26.3-27.8) 27.5 (26.7-28.2)" Great! So now our data looks much better but we need to add back our header and we would like this to only be a single line to make it easy to transform our data into a table or table-like object.
As a reminder, here is what our header used to look like:
We will add a new header after first dealing with some additional spacing and missing data issues.
First let’s try splitting our header-less data into columns based on spaces using the str_split() function, where we specify that we are splitting the data based on the pattern of a space (the space is included in quotes):
Here, we will take a look at just the first 10 lines using the [1:10,] syntax:
[,1] [,2] [,3] [,4] [,5]
[1,] "Men" "20.2" "(17.8-22.7)" "19.7" "(17.2-22.2)"
[2,] "Afghanistan" "Women" "20.6" "(18.4-22.8)" "20.1"
[3,] "Men" "25.2" "(23.9-26.5)" "25.0" "(23.7-26.4)"
[4,] "Albania" "Women" "26.0" "(24.1-27.9)" "26.1"
[5,] "Men" "22.1" "(20.8-23.3)" "21.8" "(20.5-23.1)"
[6,] "Algeria" "Women" "24.0" "(22.2-25.7)" "23.3"
[7,] "Men" "33.7" "(32.7-34.7)" "32.6" "(31.7-33.5)"
[8,] "American" "Samoa" "Women" "34.3" "(33.1-35.6)"
[9,] "Men" "25.0" "(22.5-27.4)" "25.3" "(22.8-27.7)"
[10,] "Andorra" "Women" "25.2" "(22.0-28.4)" "25.4"
[,6] [,7] [,8] [,9] [,10]
[1,] "22.4" "(20.0-25.0)" "22.8" "(20.3-25.3)" "22.5"
[2,] "(17.8-22.4)" "23.2" "(20.9-25.4)" "24.4" "(23.3-25.4)"
[3,] "25.4" "(24.0-26.7)" "27.0" "(26.0-27.9)" "26.9"
[4,] "(24.1-28.1)" "25.9" "(23.9-27.8)" "26.0" "(24.8-27.2)"
[5,] "22.3" "(21.0-23.6)" "25.1" "(24.5-25.7)" "24.8"
[6,] "(21.4-25.1)" "24.8" "(22.9-26.6)" "27.4" "(26.7-28.0)"
[7,] "34.0" "(32.9-35.1)" "34.3" "(33.0-35.6)" "34.6"
[8,] "34.1" "(33.0-35.2)" "34.4" "(33.0-35.8)" "35.3"
[9,] "25.0" "(22.5-27.4)" "26.8" "(24.4-29.2)" "26.8"
[10,] "(22.2-28.7)" "25.2" "(22.0-28.4)" "25.3" "(22.1-28.6)"
[,11] [,12] [,13] [,14] [,15]
[1,] "(20.0-25.0)" "23.6" "(21.0-26.1)" "" ""
[2,] "23.6" "(22.5-24.8)" "26.3" "(25.1-27.4)" ""
[3,] "(25.9-27.9)" "27.0" "(26.0-28.0)" "" ""
[4,] "26.2" "(24.8-27.5)" "25.9" "(24.6-27.2)" ""
[5,] "(24.1-25.4)" "25.2" "(24.6-25.9)" "" ""
[6,] "27.0" "(26.3-27.8)" "27.5" "(26.7-28.2)" ""
[7,] "(33.1-35.9)" "34.2" "(32.9-35.6)" "" ""
[8,] "(33.7-36.9)" "35.0" "(33.1-36.9)" "35.4" "(33.7-37.1)"
[9,] "(24.3-29.2)" "26.8" "(24.3-29.3)" "" ""
[10,] "25.2" "(21.9-28.5)" "25.3" "(22.1-28.6)" ""
[,16] [,17] [,18]
[1,] "" "" ""
[2,] "" "" ""
[3,] "" "" ""
[4,] "" "" ""
[5,] "" "" ""
[6,] "" "" ""
[7,] "" "" ""
[8,] "" "" ""
[9,] "" "" ""
[10,] "" "" "" This almost worked, but unfortunately country names that have spaces will be a problem. For example, we can see that “American Samoa” has been divided into two columns and all subsequent columns are shifted.
Let’s try this another way.
As we can see by looking at our rural_ubran data, the country information is only present when the sex is Women. Let’s try to extract the country information from the rows that contain data about Women.
[1] "Men 20.2 (17.8-22.7) 19.7 (17.2-22.2) 22.4 (20.0-25.0) 22.8 (20.3-25.3) 22.5 (20.0-25.0) 23.6 (21.0-26.1)"
[2] "Afghanistan Women 20.6 (18.4-22.8) 20.1 (17.8-22.4) 23.2 (20.9-25.4) 24.4 (23.3-25.4) 23.6 (22.5-24.8) 26.3 (25.1-27.4)"
[3] "Men 25.2 (23.9-26.5) 25.0 (23.7-26.4) 25.4 (24.0-26.7) 27.0 (26.0-27.9) 26.9 (25.9-27.9) 27.0 (26.0-28.0)"
[4] "Albania Women 26.0 (24.1-27.9) 26.1 (24.1-28.1) 25.9 (23.9-27.8) 26.0 (24.8-27.2) 26.2 (24.8-27.5) 25.9 (24.6-27.2)"
[5] "Men 22.1 (20.8-23.3) 21.8 (20.5-23.1) 22.3 (21.0-23.6) 25.1 (24.5-25.7) 24.8 (24.1-25.4) 25.2 (24.6-25.9)"
[6] "Algeria Women 24.0 (22.2-25.7) 23.3 (21.4-25.1) 24.8 (22.9-26.6) 27.4 (26.7-28.0) 27.0 (26.3-27.8) 27.5 (26.7-28.2)" To do this, we start by noticing that the sex always starts with either a capital “W” if the sex is female. Note, we need to use a space before the “W” otherwise we will split some of the country names if the names starts with “W”.
Here, we will also introduce the concept of piping, which uses the %>% from the magrittr package, which allows us to pipe the output from one function to the input for another function. This is really useful when we have multiple steps, which we will show soon.
First, we will select just the data for Women and Men separately:
Women <- str_subset(string = rural_urban,
pattern = "Women")
Men <- str_subset(string = rural_urban,
pattern = "Men")
head(Women)[1] "Afghanistan Women 20.6 (18.4-22.8) 20.1 (17.8-22.4) 23.2 (20.9-25.4) 24.4 (23.3-25.4) 23.6 (22.5-24.8) 26.3 (25.1-27.4)"
[2] "Albania Women 26.0 (24.1-27.9) 26.1 (24.1-28.1) 25.9 (23.9-27.8) 26.0 (24.8-27.2) 26.2 (24.8-27.5) 25.9 (24.6-27.2)"
[3] "Algeria Women 24.0 (22.2-25.7) 23.3 (21.4-25.1) 24.8 (22.9-26.6) 27.4 (26.7-28.0) 27.0 (26.3-27.8) 27.5 (26.7-28.2)"
[4] "American Samoa Women 34.3 (33.1-35.6) 34.1 (33.0-35.2) 34.4 (33.0-35.8) 35.3 (33.7-36.9) 35.0 (33.1-36.9) 35.4 (33.7-37.1)"
[5] "Andorra Women 25.2 (22.0-28.4) 25.4 (22.2-28.7) 25.2 (22.0-28.4) 25.3 (22.1-28.6) 25.2 (21.9-28.5) 25.3 (22.1-28.6)"
[6] "Angola Women 21.3 (18.0-24.6) 20.9 (17.6-24.3) 22.7 (19.3-26.0) 24.4 (21.2-27.7) 23.3 (20.0-26.7) 25.7 (22.4-29.1)" [1] "Men 20.2 (17.8-22.7) 19.7 (17.2-22.2) 22.4 (20.0-25.0) 22.8 (20.3-25.3) 22.5 (20.0-25.0) 23.6 (21.0-26.1)"
[2] "Men 25.2 (23.9-26.5) 25.0 (23.7-26.4) 25.4 (24.0-26.7) 27.0 (26.0-27.9) 26.9 (25.9-27.9) 27.0 (26.0-28.0)"
[3] "Men 22.1 (20.8-23.3) 21.8 (20.5-23.1) 22.3 (21.0-23.6) 25.1 (24.5-25.7) 24.8 (24.1-25.4) 25.2 (24.6-25.9)"
[4] "Men 33.7 (32.7-34.7) 32.6 (31.7-33.5) 34.0 (32.9-35.1) 34.3 (33.0-35.6) 34.6 (33.1-35.9) 34.2 (32.9-35.6)"
[5] "Men 25.0 (22.5-27.4) 25.3 (22.8-27.7) 25.0 (22.5-27.4) 26.8 (24.4-29.2) 26.8 (24.3-29.2) 26.8 (24.3-29.3)"
[6] "Men 20.5 (17.9-23.1) 20.2 (17.6-22.9) 21.4 (18.8-24.0) 22.6 (20.0-25.1) 22.0 (19.4-24.6) 23.2 (20.6-25.9)"And then split the Women data based on the pattern pattern = " Women" (Note that the space before the word “Women” is included - so that it is not within the country data). Here we take our Women dataset and pipe it into the first argument of the str_split() function, where then we can split on the pattern " Women". We then pipe through the unlist() function so that the results is a vector instead of a list. We assign this new result to the country_split object.
[1] "Afghanistan"
[2] " 20.6 (18.4-22.8) 20.1 (17.8-22.4) 23.2 (20.9-25.4) 24.4 (23.3-25.4) 23.6 (22.5-24.8) 26.3 (25.1-27.4)"
[3] "Albania"
[4] " 26.0 (24.1-27.9) 26.1 (24.1-28.1) 25.9 (23.9-27.8) 26.0 (24.8-27.2) 26.2 (24.8-27.5) 25.9 (24.6-27.2)"
[5] "Algeria"
[6] " 24.0 (22.2-25.7) 23.3 (21.4-25.1) 24.8 (22.9-26.6) 27.4 (26.7-28.0) 27.0 (26.3-27.8) 27.5 (26.7-28.2)"However, we can also do the same thing in a single step, where we replace the Women dataset with the code we used to create that dataset. So first, we select the “Women” and then we split the data by " Women" to create a new dataset all within the same chunk of code.
country_split <-
str_subset(string = rural_urban,
pattern = "Women") %>%
stringr::str_split(pattern= " Women") %>%
unlist()
head(country_split)[1] "Afghanistan"
[2] " 20.6 (18.4-22.8) 20.1 (17.8-22.4) 23.2 (20.9-25.4) 24.4 (23.3-25.4) 23.6 (22.5-24.8) 26.3 (25.1-27.4)"
[3] "Albania"
[4] " 26.0 (24.1-27.9) 26.1 (24.1-28.1) 25.9 (23.9-27.8) 26.0 (24.8-27.2) 26.2 (24.8-27.5) 25.9 (24.6-27.2)"
[5] "Algeria"
[6] " 24.0 (22.2-25.7) 23.3 (21.4-25.1) 24.8 (22.9-26.6) 27.4 (26.7-28.0) 27.0 (26.3-27.8) 27.5 (26.7-28.2)"Now, we can see that Country is always the odd rows and the BMI data for Women is the even rows. We can select the odd rows using the Modulus (Remainder from division) operator. This operator looks like this %%. All odd values such as 3, 7, or 15 when divided by 2 would leave a remainder of 1. We will select just these rows using the filter() function of dplyr.
We can take a look to make sure that all the country names look as expected. You can scroll through the country names if you hover over the data.
country_split
1 Afghanistan
2 Albania
3 Algeria
4 American Samoa
5 Andorra
6 Angola
7 Antigua and Barbuda
8 Argentina
9 Armenia
10 Australia
11 Austria
12 Azerbaijan
13 Bahamas
14 Bahrain
15 Bangladesh
16 Barbados
17 Belarus
18 Belgium
19 Belize
20 Benin
21 Bermuda
22 Bhutan
23 Bolivia
24 Bosnia and Herzegovina
25 Botswana
26 Brazil
27 Brunei Darussalam
28 Bulgaria
29 Burkina Faso
30 Burundi
31 Cabo Verde
32 Cambodia
33 Cameroon
34 Canada
35 Central African Republic
36 Chad
37 Chile
38 China
39 China (Hong Kong SAR)
40 Colombia
41 Comoros
42 Congo
43 Cook Islands
44 Costa Rica
45 Cote d'Ivoire
46 Croatia
47 Cuba
48 Cyprus
49 Czech Republic
50 Denmark
51 Djibouti
52 Dominica
53 Dominican Republic
54 DR Congo
55 Ecuador
56 Egypt
57 El Salvador
58 Equatorial Guinea
59 Eritrea
60 Estonia
61 Ethiopia
62 Fiji
63 Finland
64 France
65 French Polynesia
66 Gabon
67 Gambia
68 Georgia
69 Germany
70 Ghana
71 Greece
72 Greenland
73 Grenada
74 Guatemala
75 Guinea
76 Guinea Bissau
77 Guyana
78 Haiti
79 Honduras
80 Hungary
81 Iceland
82 India
83 Indonesia
84 Iran
85 Iraq
86 Ireland
87 Israel
88 Italy
89 Jamaica
90 Japan
91 Jordan
92 Kazakhstan
93 Kenya
94 Kiribati
95 Kuwait
96 Kyrgyzstan
97 Lao PDR
98 Latvia
99 Lebanon
100 Lesotho
101 Liberia
102 Libya
103 Lithuania
104 Luxembourg
105 Macedonia (TFYR)
106 Madagascar
107 Malawi
108 Malaysia
109 Maldives
110 Mali
111 Malta
112 Marshall Islands
113 Mauritania
114 Mauritius
115 Mexico
116 Micronesia (Federated States of)
117 Moldova
118 Mongolia
119 Montenegro
120 Morocco
121 Mozambique
122 Myanmar
123 Namibia
124 Nauru
125 Nepal
126 Netherlands
127 New Zealand
128 Nicaragua
129 Niger
130 Nigeria
131 Niue
132 North Korea
133 Norway
134 Occupied Palestinian Territory
135 Oman
136 Pakistan
137 Palau
138 Panama
139 Papua New Guinea
140 Paraguay
141 Peru
142 Philippines
143 Poland
144 Portugal
145 Puerto Rico
146 Qatar
147 Romania
148 Russian Federation
149 Rwanda
150 Saint Kitts and Nevis
151 Saint Lucia
152 Saint Vincent and the Grenadines
153 Samoa
154 Sao Tome and Principe
155 Saudi Arabia
156 Senegal
157 Serbia
158 Seychelles
159 Sierra Leone
160 Singapore
161 Slovakia
162 Slovenia
163 Solomon Islands
164 Somalia
165 South Africa
166 South Korea
167 Spain
168 Sri Lanka
169 Sudan (former)
170 Suriname
171 Swaziland
172 Sweden
173 Switzerland
174 Syrian Arab Republic
175 Taiwan
176 Tajikistan
177 Tanzania
178 Thailand
179 Timor-Leste
180 Togo
181 Tokelau
182 Tonga
183 Trinidad and Tobago
184 Tunisia
185 Turkey
186 Turkmenistan
187 Tuvalu
188 Uganda
189 Ukraine
190 United Arab Emirates
191 United Kingdom
192 United States of America
193 Uruguay
194 Uzbekistan
195 Vanuatu
196 Venezuela
197 Viet Nam
198 Yemen
199 Zambia
200 ZimbabweLooks good! We have 200 countries represented in our dataset.
Great! Now we have a list of the countries that can be used for both the male and female data.
Remember the even rows are the rows with the BMI data for women. These row values have a remainder of 0 when divided by 2. Let’s select that data too.
# A tibble: 6 x 1
country_split
<chr>
1 " 20.6 (18.4-22.8) 20.1 (17.8-22.4) 23.2 (20.9-25.4) 24.4 (23.3-25.4) 23.6 (2…
2 " 26.0 (24.1-27.9) 26.1 (24.1-28.1) 25.9 (23.9-27.8) 26.0 (24.8-27.2) 26.2 (2…
3 " 24.0 (22.2-25.7) 23.3 (21.4-25.1) 24.8 (22.9-26.6) 27.4 (26.7-28.0) 27.0 (2…
4 " 34.3 (33.1-35.6) 34.1 (33.0-35.2) 34.4 (33.0-35.8) 35.3 (33.7-36.9) 35.0 (3…
5 " 25.2 (22.0-28.4) 25.4 (22.2-28.7) 25.2 (22.0-28.4) 25.3 (22.1-28.6) 25.2 (2…
6 " 21.3 (18.0-24.6) 20.9 (17.6-24.3) 22.7 (19.3-26.0) 24.4 (21.2-27.7) 23.3 (2…This looks similar to our Men data (although slightly different):
[1] "Men 20.2 (17.8-22.7) 19.7 (17.2-22.2) 22.4 (20.0-25.0) 22.8 (20.3-25.3) 22.5 (20.0-25.0) 23.6 (21.0-26.1)"
[2] "Men 25.2 (23.9-26.5) 25.0 (23.7-26.4) 25.4 (24.0-26.7) 27.0 (26.0-27.9) 26.9 (25.9-27.9) 27.0 (26.0-28.0)"
[3] "Men 22.1 (20.8-23.3) 21.8 (20.5-23.1) 22.3 (21.0-23.6) 25.1 (24.5-25.7) 24.8 (24.1-25.4) 25.2 (24.6-25.9)"
[4] "Men 33.7 (32.7-34.7) 32.6 (31.7-33.5) 34.0 (32.9-35.1) 34.3 (33.0-35.6) 34.6 (33.1-35.9) 34.2 (32.9-35.6)"
[5] "Men 25.0 (22.5-27.4) 25.3 (22.8-27.7) 25.0 (22.5-27.4) 26.8 (24.4-29.2) 26.8 (24.3-29.2) 26.8 (24.3-29.3)"
[6] "Men 20.5 (17.9-23.1) 20.2 (17.6-22.9) 21.4 (18.8-24.0) 22.6 (20.0-25.1) 22.0 (19.4-24.6) 23.2 (20.6-25.9)"Notice how the data in each row is in quotes. This is because these are character strings. We can check this by using the base class() function.
[1] "character"Whereas the Women_BMI data is a tibble.
[1] "tbl_df" "tbl" "data.frame"Let’s make the Men data into a tibble as well:
It’s always a good idea to check that your data objects are the size you expect when wrangling. We can do so with the dim() function, which shows us the dimensions of data objects.
[1] 200 1[1] 200 1# A tibble: 6 x 1
country_split
<chr>
1 " 20.6 (18.4-22.8) 20.1 (17.8-22.4) 23.2 (20.9-25.4) 24.4 (23.3-25.4) 23.6 (2…
2 " 26.0 (24.1-27.9) 26.1 (24.1-28.1) 25.9 (23.9-27.8) 26.0 (24.8-27.2) 26.2 (2…
3 " 24.0 (22.2-25.7) 23.3 (21.4-25.1) 24.8 (22.9-26.6) 27.4 (26.7-28.0) 27.0 (2…
4 " 34.3 (33.1-35.6) 34.1 (33.0-35.2) 34.4 (33.0-35.8) 35.3 (33.7-36.9) 35.0 (3…
5 " 25.2 (22.0-28.4) 25.4 (22.2-28.7) 25.2 (22.0-28.4) 25.3 (22.1-28.6) 25.2 (2…
6 " 21.3 (18.0-24.6) 20.9 (17.6-24.3) 22.7 (19.3-26.0) 24.4 (21.2-27.7) 23.3 (2…# A tibble: 6 x 1
Men
<chr>
1 Men 20.2 (17.8-22.7) 19.7 (17.2-22.2) 22.4 (20.0-25.0) 22.8 (20.3-25.3) 22.5 …
2 Men 25.2 (23.9-26.5) 25.0 (23.7-26.4) 25.4 (24.0-26.7) 27.0 (26.0-27.9) 26.9 …
3 Men 22.1 (20.8-23.3) 21.8 (20.5-23.1) 22.3 (21.0-23.6) 25.1 (24.5-25.7) 24.8 …
4 Men 33.7 (32.7-34.7) 32.6 (31.7-33.5) 34.0 (32.9-35.1) 34.3 (33.0-35.6) 34.6 …
5 Men 25.0 (22.5-27.4) 25.3 (22.8-27.7) 25.0 (22.5-27.4) 26.8 (24.4-29.2) 26.8 …
6 Men 20.5 (17.9-23.1) 20.2 (17.6-22.9) 21.4 (18.8-24.0) 22.6 (20.0-25.1) 22.0 …Great! There are 200 rows of character strings with the BMI data like we expected to match the 200 countries in our data.
We do however need to add the word " Women" back to the Women_BMI data. To do this we will use the glue() function of the glue package.
Here is a simple example of how this function glues two things together.
the cow jumped over the moonNow let’s add the word “Women” to each row of the Women_BMI data using glue.
First let’s pull the country_split data out as a vector to make this simpler when we glue. To do this we will use the pull() function of the dplyr package. This function will allow us to extract or pull out the single variable called country_split which contains the BMI values for women from within the tibble called Women_BMI which we just created by grabbing only the even rows of our original country_split tibble.
Then, let’s use the first() function of the dplyr package to get a reminder of what just the first row of the country_split column of the Women_BMI data looks like.
[1] " 20.6 (18.4-22.8) 20.1 (17.8-22.4) 23.2 (20.9-25.4) 24.4 (23.3-25.4) 23.6 (22.5-24.8) 26.3 (25.1-27.4)"Now we will paste the word "Women" in the beginning of each value in the vector using the glue() function.
Note that there will be a space in between. This is because the space is already in front of the first BMI value.
Women 20.6 (18.4-22.8) 20.1 (17.8-22.4) 23.2 (20.9-25.4) 24.4 (23.3-25.4) 23.6 (22.5-24.8) 26.3 (25.1-27.4)
Women 26.0 (24.1-27.9) 26.1 (24.1-28.1) 25.9 (23.9-27.8) 26.0 (24.8-27.2) 26.2 (24.8-27.5) 25.9 (24.6-27.2)
Women 24.0 (22.2-25.7) 23.3 (21.4-25.1) 24.8 (22.9-26.6) 27.4 (26.7-28.0) 27.0 (26.3-27.8) 27.5 (26.7-28.2)
Women 34.3 (33.1-35.6) 34.1 (33.0-35.2) 34.4 (33.0-35.8) 35.3 (33.7-36.9) 35.0 (33.1-36.9) 35.4 (33.7-37.1)
Women 25.2 (22.0-28.4) 25.4 (22.2-28.7) 25.2 (22.0-28.4) 25.3 (22.1-28.6) 25.2 (21.9-28.5) 25.3 (22.1-28.6)
Women 21.3 (18.0-24.6) 20.9 (17.6-24.3) 22.7 (19.3-26.0) 24.4 (21.2-27.7) 23.3 (20.0-26.7) 25.7 (22.4-29.1)You could also do all of this in one step but its less obvious what is happening.
Here we can pull the country_split data out of the Women_BMI object and then paste "Women" in front of each string in one command.
If we try splitting our data by space again, will it have the expected number of columns? What about the rows that contain na* values?
Let’s just take the Men data that contains na* values. This column is called rural_urban.
# A tibble: 5 x 1
Men
<chr>
1 Men 24.3 (21.8-26.8) na* 24.3 (21.8-26.8) 26.7 (24.3-29.2) na* 26.7 (24.3-29.…
2 Men 22.4 (21.4-23.3) 21.5 (20.2-22.8) 22.4 (21.5-23.4) 24.8 (23.6-26.1) na* 2…
3 Men 32.6 (32.0-33.2) na* 32.6 (32.0-33.2) 32.9 (31.7-34.1) na* 32.9 (31.7-34.…
4 Men 22.8 (22.2-23.3) na* 22.8 (22.2-23.3) 24.4 (23.6-25.1) na* 24.4 (23.6-25.…
5 Men 29.6 (28.1-31.1) 29.6 (28.1-31.1) na* 32.3 (31.4-33.2) 32.3 (31.4-33.2) n…Now we can try splitting by a space
# A tibble: 5 x 12
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12
<chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 Men 24.3 (21.8… na* 24.3 (21.8… 26.7 (24.3… na* 26.7 (24.… ""
2 Men 22.4 (21.4… 21.5 (20.2… 22.4 (21.5… 24.8 (23.6… na* 24.8 "(23.…
3 Men 32.6 (32.0… na* 32.6 (32.0… 32.9 (31.7… na* 32.9 (31.… ""
4 Men 22.8 (22.2… na* 22.8 (22.2… 24.4 (23.6… na* 24.4 (23.… ""
5 Men 29.6 (28.1… 29.6 (28.1… na* 32.3 (31.4… 32.3 (31.… na* "" So close! Notice that the "na*" values have shifted the subsequent values within the columns because typically there is a space between the BMI and the credible intervals. Here we can see this data in our original pdf:
We need to replace our na* values with something that includes a space so that when we separate our data by space we will have two values instead of one when we have an na*. Therefore, na* na* should work.
So, in order to use the functions within the stringr package, our Men data needs to be of character class, so let’s check that now.
[1] "character"The as.character() function allows us to change factor data within to the character class. We could use this if we needed it, but it looks like our data is already character.
We are also changing the structure so that the data is not within a column called Men.
# A tibble: 6 x 1
Men
<chr>
1 Men 20.2 (17.8-22.7) 19.7 (17.2-22.2) 22.4 (20.0-25.0) 22.8 (20.3-25.3) 22.5 …
2 Men 25.2 (23.9-26.5) 25.0 (23.7-26.4) 25.4 (24.0-26.7) 27.0 (26.0-27.9) 26.9 …
3 Men 22.1 (20.8-23.3) 21.8 (20.5-23.1) 22.3 (21.0-23.6) 25.1 (24.5-25.7) 24.8 …
4 Men 33.7 (32.7-34.7) 32.6 (31.7-33.5) 34.0 (32.9-35.1) 34.3 (33.0-35.6) 34.6 …
5 Men 25.0 (22.5-27.4) 25.3 (22.8-27.7) 25.0 (22.5-27.4) 26.8 (24.4-29.2) 26.8 …
6 Men 20.5 (17.9-23.1) 20.2 (17.6-22.9) 21.4 (18.8-24.0) 22.6 (20.0-25.1) 22.0 …[1] "Men 20.2 (17.8-22.7) 19.7 (17.2-22.2) 22.4 (20.0-25.0) 22.8 (20.3-25.3) 22.5 (20.0-25.0) 23.6 (21.0-26.1)"
[2] "Men 25.2 (23.9-26.5) 25.0 (23.7-26.4) 25.4 (24.0-26.7) 27.0 (26.0-27.9) 26.9 (25.9-27.9) 27.0 (26.0-28.0)"
[3] "Men 22.1 (20.8-23.3) 21.8 (20.5-23.1) 22.3 (21.0-23.6) 25.1 (24.5-25.7) 24.8 (24.1-25.4) 25.2 (24.6-25.9)"
[4] "Men 33.7 (32.7-34.7) 32.6 (31.7-33.5) 34.0 (32.9-35.1) 34.3 (33.0-35.6) 34.6 (33.1-35.9) 34.2 (32.9-35.6)"
[5] "Men 25.0 (22.5-27.4) 25.3 (22.8-27.7) 25.0 (22.5-27.4) 26.8 (24.4-29.2) 26.8 (24.3-29.2) 26.8 (24.3-29.3)"
[6] "Men 20.5 (17.9-23.1) 20.2 (17.6-22.9) 21.4 (18.8-24.0) 22.6 (20.0-25.1) 22.0 (19.4-24.6) 23.2 (20.6-25.9)"Men <-
stringr::str_replace_all(string = Men,
pattern = "na\\*",
replacement = "na\\* na\\*")
# The \\ are necessary because the * is a special character
# The * would typically indicate any possible value,
# but here we actually want a "*" instead
# Thus the double backslash does that for us
# Here we are replacing all occurences of the na* values
#(thus str_replace_all instead of str_replace) with na* na*. Let’s check that it worked…
[1] "Men 20.7 (19.4-22.1) 20.3 (18.9-21.7) 21.7 (20.3-23.1) 22.7 (22.0-23.3) 22.1 (21.3-22.8) 23.4 (22.7-24.1)"
[2] "Men 24.3 (21.8-26.8) na* na* 24.3 (21.8-26.8) 26.7 (24.3-29.2) na* na* 26.7 (24.3-29.2)"
[3] "Men 20.6 (19.1-22.1) 20.3 (18.7-21.8) 23.1 (21.5-24.6) 23.5 (22.8-24.3) 23.1 (22.2-23.9) 24.2 (23.3-25.2)"
[4] "Men 23.3 (21.0-25.5) 22.9 (20.6-25.2) 23.6 (21.3-25.9) 26.1 (23.9-28.4) 25.3 (23.1-27.5) 26.5 (24.2-28.8)"
[5] "Men 24.9 (23.5-26.2) 24.8 (23.4-26.2) 25.0 (23.6-26.4) 26.8 (25.7-27.8) 26.7 (25.7-27.8) 26.8 (25.7-27.9)"
[6] "Men 20.6 (19.2-22.0) 20.3 (18.8-21.7) 21.4 (20.0-22.9) 22.6 (21.9-23.4) 21.9 (20.9-22.8) 23.2 (22.3-24.0)"
[7] "Men 23.3 (22.5-24.0) 22.2 (21.4-23.0) 23.7 (22.9-24.5) 26.2 (25.7-26.7) 25.0 (24.3-25.6) 26.4 (25.9-26.9)"
[8] "Men 24.0 (22.5-25.4) 23.7 (22.1-25.1) 24.2 (22.7-25.7) 27.1 (26.4-27.7) 26.6 (26.0-27.3) 27.2 (26.5-27.9)"
[9] "Men 25.0 (23.8-26.4) 24.9 (23.5-26.4) 25.1 (23.8-26.5) 27.0 (25.8-28.1) 27.0 (25.7-28.3) 27.0 (25.7-28.2)"
[10] "Men 20.2 (18.8-21.5) 20.0 (18.6-21.4) 21.4 (19.9-22.9) 22.2 (21.4-23.0) 21.8 (21.0-22.6) 23.1 (22.1-24.1)"
[11] "Men 20.2 (17.7-22.6) 20.1 (17.6-22.5) 21.2 (18.7-23.6) 22.1 (19.8-24.5) 22.0 (19.6-24.3) 23.2 (20.8-25.6)"Great!
Now for the Women data object
[1] "glue" "character"# the female data is already character type
Women <-
stringr::str_replace_all(string = Women,
pattern = "na\\*",
replacement = "na\\* na\\*")
Women[20:30] [1] "Women 20.6 (19.1-22.0) 20.1 (18.5-21.6) 21.8 (20.3-23.3) 24.3 (23.5-25.0) 23.2 (22.5-24.0) 25.5 (24.7-26.3)"
[2] "Women 25.4 (22.2-28.6) na* na* 25.4 (22.2-28.6) 28.5 (25.3-31.6) na* na* 28.5 (25.3-31.6)"
[3] "Women 20.7 (18.4-22.9) 20.3 (18.0-22.6) 23.0 (20.8-25.3) 24.6 (23.7-25.5) 23.8 (22.7-24.9) 25.9 (24.7-27.0)"
[4] "Women 23.8 (22.3-25.3) 23.0 (21.4-24.5) 24.6 (23.1-26.1) 27.9 (26.5-29.3) 27.0 (25.6-28.5) 28.3 (26.8-29.7)"
[5] "Women 25.5 (23.5-27.4) 25.7 (23.7-27.7) 25.1 (23.1-27.0) 25.7 (24.4-27.1) 26.0 (24.6-27.5) 25.3 (23.9-26.8)"
[6] "Women 24.2 (22.2-26.2) 23.7 (21.5-25.8) 25.6 (23.5-27.7) 26.2 (25.3-27.1) 25.0 (23.8-26.2) 27.0 (26.0-28.2)"
[7] "Women 24.0 (23.2-25.0) 23.1 (22.1-24.1) 24.4 (23.5-25.4) 26.9 (26.2-27.5) 26.4 (25.7-27.2) 26.9 (26.3-27.6)"
[8] "Women 23.6 (21.4-25.6) 22.9 (20.7-25.0) 24.0 (21.8-26.1) 27.4 (26.6-28.1) 27.0 (26.1-27.8) 27.5 (26.7-28.2)"
[9] "Women 25.2 (23.4-27.1) 25.5 (23.4-27.6) 25.0 (23.1-26.9) 25.6 (24.1-27.1) 26.0 (24.2-27.8) 25.5 (23.9-27.1)"
[10] "Women 19.9 (18.7-21.2) 19.6 (18.4-20.9) 22.0 (20.7-23.3) 22.2 (21.4-23.1) 21.1 (20.2-22.0) 24.7 (23.7-25.7)"
[11] "Women 19.6 (17.5-21.7) 19.5 (17.4-21.6) 21.7 (19.5-23.8) 21.1 (20.3-22.0) 20.7 (19.8-21.6) 24.0 (23.0-24.9)"Great, now we can split our data by spaces.
We will split our data finally by spaces using the str_split() function of the stringr package. We can check the rows with NA values by looking at rows 20-30.
# A tibble: 11 x 13
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13
<chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 Men 20.7 (19.… 20.3 (18.… 21.7 (20.… 22.7 (22.… 22.1 (21.… 23.4 (22.…
2 Men 24.3 (21.… na* na* 24.3 (21.… 26.7 (24.… na* na* 26.7 (24.…
3 Men 20.6 (19.… 20.3 (18.… 23.1 (21.… 23.5 (22.… 23.1 (22.… 24.2 (23.…
4 Men 23.3 (21.… 22.9 (20.… 23.6 (21.… 26.1 (23.… 25.3 (23.… 26.5 (24.…
5 Men 24.9 (23.… 24.8 (23.… 25.0 (23.… 26.8 (25.… 26.7 (25.… 26.8 (25.…
6 Men 20.6 (19.… 20.3 (18.… 21.4 (20.… 22.6 (21.… 21.9 (20.… 23.2 (22.…
7 Men 23.3 (22.… 22.2 (21.… 23.7 (22.… 26.2 (25.… 25.0 (24.… 26.4 (25.…
8 Men 24.0 (22.… 23.7 (22.… 24.2 (22.… 27.1 (26.… 26.6 (26.… 27.2 (26.…
9 Men 25.0 (23.… 24.9 (23.… 25.1 (23.… 27.0 (25.… 27.0 (25.… 27.0 (25.…
10 Men 20.2 (18.… 20.0 (18.… 21.4 (19.… 22.2 (21.… 21.8 (21.… 23.1 (22.…
11 Men 20.2 (17.… 20.1 (17.… 21.2 (18.… 22.1 (19.… 22.0 (19.… 23.2 (20.…Great! It looks good!
Now for the Women data.
# A tibble: 6 x 13
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13
<chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 Women 20.6 (18.4… 20.1 (17.… 23.2 (20.… 24.4 (23.… 23.6 (22.… 26.3 (25.…
2 Women 26.0 (24.1… 26.1 (24.… 25.9 (23.… 26.0 (24.… 26.2 (24.… 25.9 (24.…
3 Women 24.0 (22.2… 23.3 (21.… 24.8 (22.… 27.4 (26.… 27.0 (26.… 27.5 (26.…
4 Women 34.3 (33.1… 34.1 (33.… 34.4 (33.… 35.3 (33.… 35.0 (33.… 35.4 (33.…
5 Women 25.2 (22.0… 25.4 (22.… 25.2 (22.… 25.3 (22.… 25.2 (21.… 25.3 (22.…
6 Women 21.3 (18.0… 20.9 (17.… 22.7 (19.… 24.4 (21.… 23.3 (20.… 25.7 (22.…# A tibble: 11 x 13
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13
<chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 Women 20.6 (19.… 20.1 (18.… 21.8 (20.… 24.3 (23.… 23.2 (22.… 25.5 (24.…
2 Women 25.4 (22.… na* na* 25.4 (22.… 28.5 (25.… na* na* 28.5 (25.…
3 Women 20.7 (18.… 20.3 (18.… 23.0 (20.… 24.6 (23.… 23.8 (22.… 25.9 (24.…
4 Women 23.8 (22.… 23.0 (21.… 24.6 (23.… 27.9 (26.… 27.0 (25.… 28.3 (26.…
5 Women 25.5 (23.… 25.7 (23.… 25.1 (23.… 25.7 (24.… 26.0 (24.… 25.3 (23.…
6 Women 24.2 (22.… 23.7 (21.… 25.6 (23.… 26.2 (25.… 25.0 (23.… 27.0 (26.…
7 Women 24.0 (23.… 23.1 (22.… 24.4 (23.… 26.9 (26.… 26.4 (25.… 26.9 (26.…
8 Women 23.6 (21.… 22.9 (20.… 24.0 (21.… 27.4 (26.… 27.0 (26.… 27.5 (26.…
9 Women 25.2 (23.… 25.5 (23.… 25.0 (23.… 25.6 (24.… 26.0 (24.… 25.5 (23.…
10 Women 19.9 (18.… 19.6 (18.… 22.0 (20.… 22.2 (21.… 21.1 (20.… 24.7 (23.…
11 Women 19.6 (17.… 19.5 (17.… 21.7 (19.… 21.1 (20.… 20.7 (19.… 24.0 (23.…We can see that our na values look correct.
Looks good!
We can see from our pdf and our object called header what the header was like in the original pdf document. We need to add column names to our data based on the original data.
We’ll wait to add Country to the header until we add the country information to our data.
new_header <- c("Sex","National_BMI_1985",
"National_BMI_1985_CI", "Rural_BMI_1985",
"Rural_BMI_1985_CI", "Urban_BMI_1985",
"Urban_BMI_1985_CI", "National_BMI_2017",
"National_BMI_2017_CI","Rural_BMI_2017",
"Rural_BMI_2017_CI", "Urban_BMI_2017",
"Urban_BMI_2017_CI")Let’s change the names of our columns of our tibbles to this new header for our Men and Women data
# A tibble: 6 x 13
Sex National_BMI_19… National_BMI_19… Rural_BMI_1985 Rural_BMI_1985_…
<chr> <chr> <chr> <chr> <chr>
1 Women 20.6 (18.4-22.8) 20.1 (17.8-22.4)
2 Women 26.0 (24.1-27.9) 26.1 (24.1-28.1)
3 Women 24.0 (22.2-25.7) 23.3 (21.4-25.1)
4 Women 34.3 (33.1-35.6) 34.1 (33.0-35.2)
5 Women 25.2 (22.0-28.4) 25.4 (22.2-28.7)
6 Women 21.3 (18.0-24.6) 20.9 (17.6-24.3)
# … with 8 more variables: Urban_BMI_1985 <chr>, Urban_BMI_1985_CI <chr>,
# National_BMI_2017 <chr>, National_BMI_2017_CI <chr>, Rural_BMI_2017 <chr>,
# Rural_BMI_2017_CI <chr>, Urban_BMI_2017 <chr>, Urban_BMI_2017_CI <chr># A tibble: 6 x 13
Sex National_BMI_19… National_BMI_19… Rural_BMI_1985 Rural_BMI_1985_…
<chr> <chr> <chr> <chr> <chr>
1 Men 20.2 (17.8-22.7) 19.7 (17.2-22.2)
2 Men 25.2 (23.9-26.5) 25.0 (23.7-26.4)
3 Men 22.1 (20.8-23.3) 21.8 (20.5-23.1)
4 Men 33.7 (32.7-34.7) 32.6 (31.7-33.5)
5 Men 25.0 (22.5-27.4) 25.3 (22.8-27.7)
6 Men 20.5 (17.9-23.1) 20.2 (17.6-22.9)
# … with 8 more variables: Urban_BMI_1985 <chr>, Urban_BMI_1985_CI <chr>,
# National_BMI_2017 <chr>, National_BMI_2017_CI <chr>, Rural_BMI_2017 <chr>,
# Rural_BMI_2017_CI <chr>, Urban_BMI_2017 <chr>, Urban_BMI_2017_CI <chr>Now we will add our country data to both our Men and Women tibbles using the bind_cols() function. This will add the country as a new column to the Women data object on the left.
# A tibble: 6 x 14
country_split Sex National_BMI_19… National_BMI_19… Rural_BMI_1985
<chr> <chr> <chr> <chr> <chr>
1 Afghanistan Women 20.6 (18.4-22.8) 20.1
2 Albania Women 26.0 (24.1-27.9) 26.1
3 Algeria Women 24.0 (22.2-25.7) 23.3
4 American Sam… Women 34.3 (33.1-35.6) 34.1
5 Andorra Women 25.2 (22.0-28.4) 25.4
6 Angola Women 21.3 (18.0-24.6) 20.9
# … with 9 more variables: Rural_BMI_1985_CI <chr>, Urban_BMI_1985 <chr>,
# Urban_BMI_1985_CI <chr>, National_BMI_2017 <chr>,
# National_BMI_2017_CI <chr>, Rural_BMI_2017 <chr>, Rural_BMI_2017_CI <chr>,
# Urban_BMI_2017 <chr>, Urban_BMI_2017_CI <chr>Similarly, for the data about Men, this will add the country as a new column to the Men data object on the left.
# A tibble: 6 x 14
country_split Sex National_BMI_19… National_BMI_19… Rural_BMI_1985
<chr> <chr> <chr> <chr> <chr>
1 Afghanistan Men 20.2 (17.8-22.7) 19.7
2 Albania Men 25.2 (23.9-26.5) 25.0
3 Algeria Men 22.1 (20.8-23.3) 21.8
4 American Sam… Men 33.7 (32.7-34.7) 32.6
5 Andorra Men 25.0 (22.5-27.4) 25.3
6 Angola Men 20.5 (17.9-23.1) 20.2
# … with 9 more variables: Rural_BMI_1985_CI <chr>, Urban_BMI_1985 <chr>,
# Urban_BMI_1985_CI <chr>, National_BMI_2017 <chr>,
# National_BMI_2017_CI <chr>, Rural_BMI_2017 <chr>, Rural_BMI_2017_CI <chr>,
# Urban_BMI_2017 <chr>, Urban_BMI_2017_CI <chr>Here, we will change the variable name in the country dataset to country (currently it is called country_split). We will also introduce the concept of the assignment pipe. In this case our pipe operator looks like this %<>%. Using this additional pipe requires another package, called magrittr. The other simpler pipe options from this package are loaded with tidyverse (if you used library(tidyverse) which loads most tidyverse packages), but not this version.
The > portion of the pipe still behaves like a normal pipe, while the < portion of the pipe makes an assignment to whatever the <is pointing to, just like when we use the typical assignment operator <-.
# library(magrittr)
# We can't use the `%<>%` unless we load the magrittr package
# We have already done this but we include this for illustrative purposes.
# Here we will just use the traditional assignment strategy
Women <- dplyr::rename(Women, Country = country_split) Here, we have renamed the country_split variable to be called Country.
Here, we reassign Men using the pipe strategy.
We have renamed the country_split variable to to be called Country. We have also reassigned Men to the same data object. which has the Country variable renamed.
Let’s take a look at our data for Men and Women:
# A tibble: 6 x 14
Country Sex National_BMI_19… National_BMI_19… Rural_BMI_1985
<chr> <chr> <chr> <chr> <chr>
1 Afghan… Women 20.6 (18.4-22.8) 20.1
2 Albania Women 26.0 (24.1-27.9) 26.1
3 Algeria Women 24.0 (22.2-25.7) 23.3
4 Americ… Women 34.3 (33.1-35.6) 34.1
5 Andorra Women 25.2 (22.0-28.4) 25.4
6 Angola Women 21.3 (18.0-24.6) 20.9
# … with 9 more variables: Rural_BMI_1985_CI <chr>, Urban_BMI_1985 <chr>,
# Urban_BMI_1985_CI <chr>, National_BMI_2017 <chr>,
# National_BMI_2017_CI <chr>, Rural_BMI_2017 <chr>, Rural_BMI_2017_CI <chr>,
# Urban_BMI_2017 <chr>, Urban_BMI_2017_CI <chr># A tibble: 6 x 14
Country Sex National_BMI_19… National_BMI_19… Rural_BMI_1985
<chr> <chr> <chr> <chr> <chr>
1 Afghan… Men 20.2 (17.8-22.7) 19.7
2 Albania Men 25.2 (23.9-26.5) 25.0
3 Algeria Men 22.1 (20.8-23.3) 21.8
4 Americ… Men 33.7 (32.7-34.7) 32.6
5 Andorra Men 25.0 (22.5-27.4) 25.3
6 Angola Men 20.5 (17.9-23.1) 20.2
# … with 9 more variables: Rural_BMI_1985_CI <chr>, Urban_BMI_1985 <chr>,
# Urban_BMI_1985_CI <chr>, National_BMI_2017 <chr>,
# National_BMI_2017_CI <chr>, Rural_BMI_2017 <chr>, Rural_BMI_2017_CI <chr>,
# Urban_BMI_2017 <chr>, Urban_BMI_2017_CI <chr>Now that both of these tibbles have the same structure, we can combine our Men and Women data using the bind_rows() function of dplyr. This will combine all of the rows based on their similar column structure. All of the male data will be first and all the female data will be second.
Let’s check the size of our BMI data… it should have 400 rows (obs). We can see the end of the BMI data using the tail() function:
tibble [400 × 14] (S3: tbl_df/tbl/data.frame)
$ Country : chr [1:400] "Afghanistan" "Albania" "Algeria" "American Samoa" ...
$ Sex : chr [1:400] "Men" "Men" "Men" "Men" ...
$ National_BMI_1985 : chr [1:400] "20.2" "25.2" "22.1" "33.7" ...
$ National_BMI_1985_CI: chr [1:400] "(17.8-22.7)" "(23.9-26.5)" "(20.8-23.3)" "(32.7-34.7)" ...
$ Rural_BMI_1985 : chr [1:400] "19.7" "25.0" "21.8" "32.6" ...
$ Rural_BMI_1985_CI : chr [1:400] "(17.2-22.2)" "(23.7-26.4)" "(20.5-23.1)" "(31.7-33.5)" ...
$ Urban_BMI_1985 : chr [1:400] "22.4" "25.4" "22.3" "34.0" ...
$ Urban_BMI_1985_CI : chr [1:400] "(20.0-25.0)" "(24.0-26.7)" "(21.0-23.6)" "(32.9-35.1)" ...
$ National_BMI_2017 : chr [1:400] "22.8" "27.0" "25.1" "34.3" ...
$ National_BMI_2017_CI: chr [1:400] "(20.3-25.3)" "(26.0-27.9)" "(24.5-25.7)" "(33.0-35.6)" ...
$ Rural_BMI_2017 : chr [1:400] "22.5" "26.9" "24.8" "34.6" ...
$ Rural_BMI_2017_CI : chr [1:400] "(20.0-25.0)" "(25.9-27.9)" "(24.1-25.4)" "(33.1-35.9)" ...
$ Urban_BMI_2017 : chr [1:400] "23.6" "27.0" "25.2" "34.2" ...
$ Urban_BMI_2017_CI : chr [1:400] "(21.0-26.1)" "(26.0-28.0)" "(24.6-25.9)" "(32.9-35.6)" ...# A tibble: 6 x 14
Country Sex National_BMI_19… National_BMI_19… Rural_BMI_1985
<chr> <chr> <chr> <chr> <chr>
1 Vanuatu Women 24.4 (22.7-26.1) 24.2
2 Venezu… Women 24.0 (22.2-25.9) 22.6
3 Viet N… Women 18.5 (17.7-19.4) 18.2
4 Yemen Women 20.5 (18.9-22.0) 20.0
5 Zambia Women 21.4 (20.3-22.6) 20.7
6 Zimbab… Women 24.6 (23.5-25.9) 24.0
# … with 9 more variables: Rural_BMI_1985_CI <chr>, Urban_BMI_1985 <chr>,
# Urban_BMI_1985_CI <chr>, National_BMI_2017 <chr>,
# National_BMI_2017_CI <chr>, Rural_BMI_2017 <chr>, Rural_BMI_2017_CI <chr>,
# Urban_BMI_2017 <chr>, Urban_BMI_2017_CI <chr>Now, let’s sort the data alphabetically by Country using the arrange() function from the dplyr package:
# A tibble: 6 x 14
Country Sex National_BMI_19… National_BMI_19… Rural_BMI_1985
<chr> <chr> <chr> <chr> <chr>
1 Afghan… Men 20.2 (17.8-22.7) 19.7
2 Afghan… Women 20.6 (18.4-22.8) 20.1
3 Albania Men 25.2 (23.9-26.5) 25.0
4 Albania Women 26.0 (24.1-27.9) 26.1
5 Algeria Men 22.1 (20.8-23.3) 21.8
6 Algeria Women 24.0 (22.2-25.7) 23.3
# … with 9 more variables: Rural_BMI_1985_CI <chr>, Urban_BMI_1985 <chr>,
# Urban_BMI_1985_CI <chr>, National_BMI_2017 <chr>,
# National_BMI_2017_CI <chr>, Rural_BMI_2017 <chr>, Rural_BMI_2017_CI <chr>,
# Urban_BMI_2017 <chr>, Urban_BMI_2017_CI <chr>Our data is looking great!
Now, we might want to make sure that our observations for each variable look the way we want. In other words, if we want to make plots about National BMI in 1985 then we would need our values to be numeric. Looking at our BMI data using str(), we can see the type of data for each of our variables listed just after variable name and the ":" colon.
tibble [400 × 14] (S3: tbl_df/tbl/data.frame)
$ Country : chr [1:400] "Afghanistan" "Afghanistan" "Albania" "Albania" ...
$ Sex : chr [1:400] "Men" "Women" "Men" "Women" ...
$ National_BMI_1985 : chr [1:400] "20.2" "20.6" "25.2" "26.0" ...
$ National_BMI_1985_CI: chr [1:400] "(17.8-22.7)" "(18.4-22.8)" "(23.9-26.5)" "(24.1-27.9)" ...
$ Rural_BMI_1985 : chr [1:400] "19.7" "20.1" "25.0" "26.1" ...
$ Rural_BMI_1985_CI : chr [1:400] "(17.2-22.2)" "(17.8-22.4)" "(23.7-26.4)" "(24.1-28.1)" ...
$ Urban_BMI_1985 : chr [1:400] "22.4" "23.2" "25.4" "25.9" ...
$ Urban_BMI_1985_CI : chr [1:400] "(20.0-25.0)" "(20.9-25.4)" "(24.0-26.7)" "(23.9-27.8)" ...
$ National_BMI_2017 : chr [1:400] "22.8" "24.4" "27.0" "26.0" ...
$ National_BMI_2017_CI: chr [1:400] "(20.3-25.3)" "(23.3-25.4)" "(26.0-27.9)" "(24.8-27.2)" ...
$ Rural_BMI_2017 : chr [1:400] "22.5" "23.6" "26.9" "26.2" ...
$ Rural_BMI_2017_CI : chr [1:400] "(20.0-25.0)" "(22.5-24.8)" "(25.9-27.9)" "(24.8-27.5)" ...
$ Urban_BMI_2017 : chr [1:400] "23.6" "26.3" "27.0" "25.9" ...
$ Urban_BMI_2017_CI : chr [1:400] "(21.0-26.1)" "(25.1-27.4)" "(26.0-28.0)" "(24.6-27.2)" ...Looks like none of our BMI data is actually numeric (or num), but of the class character (or chr). Let’s change that now.
To start, we could change these values to be numeric with 6 lines of code like this, where the double brackets[[]] and quotes are used to select the specific variable within the BMI tibble. In this case we need to use the brackets instead of pull() because we can’t start a command with pull() to reassign values in a variable.
#pull(BMI, National_BMI_1985)%<>% as.numeric #This will not work
BMI[["National_BMI_1985"]] %<>% as.numeric
BMI[["Rural_BMI_1985"]] %<>% as.numeric
BMI[["Urban_BMI_1985"]] %<>% as.numeric
BMI[["National_BMI_2017"]] %<>% as.numeric
BMI[["Rural_BMI_2017"]] %<>% as.numeric
BMI[["Urban_BMI_2017"]] %<>% as.numericAnd if we did this we would see these variables show as "num" which stands for a numeric type of variable.
tibble [400 × 14] (S3: tbl_df/tbl/data.frame)
$ Country : chr [1:400] "Afghanistan" "Afghanistan" "Albania" "Albania" ...
$ Sex : chr [1:400] "Men" "Women" "Men" "Women" ...
$ National_BMI_1985 : num [1:400] 20.2 20.6 25.2 26 22.1 24 33.7 34.3 25 25.2 ...
$ National_BMI_1985_CI: chr [1:400] "(17.8-22.7)" "(18.4-22.8)" "(23.9-26.5)" "(24.1-27.9)" ...
$ Rural_BMI_1985 : num [1:400] 19.7 20.1 25 26.1 21.8 23.3 32.6 34.1 25.3 25.4 ...
$ Rural_BMI_1985_CI : chr [1:400] "(17.2-22.2)" "(17.8-22.4)" "(23.7-26.4)" "(24.1-28.1)" ...
$ Urban_BMI_1985 : num [1:400] 22.4 23.2 25.4 25.9 22.3 24.8 34 34.4 25 25.2 ...
$ Urban_BMI_1985_CI : chr [1:400] "(20.0-25.0)" "(20.9-25.4)" "(24.0-26.7)" "(23.9-27.8)" ...
$ National_BMI_2017 : num [1:400] 22.8 24.4 27 26 25.1 27.4 34.3 35.3 26.8 25.3 ...
$ National_BMI_2017_CI: chr [1:400] "(20.3-25.3)" "(23.3-25.4)" "(26.0-27.9)" "(24.8-27.2)" ...
$ Rural_BMI_2017 : num [1:400] 22.5 23.6 26.9 26.2 24.8 27 34.6 35 26.8 25.2 ...
$ Rural_BMI_2017_CI : chr [1:400] "(20.0-25.0)" "(22.5-24.8)" "(25.9-27.9)" "(24.8-27.5)" ...
$ Urban_BMI_2017 : num [1:400] 23.6 26.3 27 25.9 25.2 27.5 34.2 35.4 26.8 25.3 ...
$ Urban_BMI_2017_CI : chr [1:400] "(21.0-26.1)" "(25.1-27.4)" "(26.0-28.0)" "(24.6-27.2)" ...Or we can use a more automated way with the map() function of the purrr package:
The map() function of the purrr package allows us to apply the as.numeric() function to all the selected columns of BMI dataset. Above, we have selected those that don’t contain CI, Sex, or Country in the column name. If you are familiar with apply(), this function is very similar.
Now we need to add the Country and Sex variables back to our new numeric dataset.
# A tibble: 6 x 8
National_BMI_19… Rural_BMI_1985 Urban_BMI_1985 National_BMI_20… Rural_BMI_2017
<dbl> <dbl> <dbl> <dbl> <dbl>
1 20.2 19.7 22.4 22.8 22.5
2 20.6 20.1 23.2 24.4 23.6
3 25.2 25 25.4 27 26.9
4 26 26.1 25.9 26 26.2
5 22.1 21.8 22.3 25.1 24.8
6 24 23.3 24.8 27.4 27
# … with 3 more variables: Urban_BMI_2017 <dbl>, Country <chr>, Sex <chr>In this case our numeric data is of a class dbl. This means double-precision floating point numbers (with decimals). The alternative numeric option is the integer class.
It is generally useful to get the data in what is called long format for other analyses, and particularly for plotting.
For a more detailed description about this, please see this case study.
Basically, the long format has more rows and fewer columns. Thus, we can combine some columns together. In our case, we can put all the different BMI data together in one long column. We will create a new column that tells us what each row of the new BMI column represents. In other words, it will tell us what the original column was.
To do this we will use the pivot_longer() function of the tidyr package. The data argument specifies what data we want to work with, the cols argument specifies the names of the columns or variables that we want to combine and change to one long format variable, and the names_to argument specifies the name for the new variable that we are creating. The : is used to indicate that we want to combine the columns or variables between National_BMI_1985 and Urban_BMI_2017.
BMI_long <-
tidyr::pivot_longer(data = BMI_numeric,
cols = National_BMI_1985:Urban_BMI_2017,
names_to = "class_BMI")
head(BMI_long)# A tibble: 6 x 4
Country Sex class_BMI value
<chr> <chr> <chr> <dbl>
1 Afghanistan Men National_BMI_1985 20.2
2 Afghanistan Men Rural_BMI_1985 19.7
3 Afghanistan Men Urban_BMI_1985 22.4
4 Afghanistan Men National_BMI_2017 22.8
5 Afghanistan Men Rural_BMI_2017 22.5
6 Afghanistan Men Urban_BMI_2017 23.6We really want to promote the idea of working with tidy data. Tidy data is data that is structured in a way that makes it especially usable. One of the rules of tidy data, is that each column or variable should have only one type of data. Currently we have two types of data in the BMI_class variable: the type of BMI measurement ( Rural, Urban, or National), and the year (1985 and 2017). We can separate these two data types of the class_BMI column using the separate() function of the tidyr package to replace the class_BMI column with two new columns.
Here, we separate the class_BMI column of the BMI_long tibble. This is based on the underscores and create two new columns. The first column will be called Region. It will contain the 1st part of the class_BMI data before the 1st underscore. The 2nd column will be called Year. It will contain the 3rd part of the data based on underscores. The middle part of the column will not be used to create any new columns.
# A tibble: 6 x 5
Country Sex Region Year value
<chr> <chr> <chr> <chr> <dbl>
1 Afghanistan Men National 1985 20.2
2 Afghanistan Men Rural 1985 19.7
3 Afghanistan Men Urban 1985 22.4
4 Afghanistan Men National 2017 22.8
5 Afghanistan Men Rural 2017 22.5
6 Afghanistan Men Urban 2017 23.6Let’s see how the dimensions of the BMI data have changed in long format:
[1] 400 8[1] 2400 5Let’s rename the value column to something more informative. We can use the rename() function of the dplyr package to do this. Notice that the value we want to replace is listed second after the = sign.
# A tibble: 6 x 5
Country Sex Region Year BMI
<chr> <chr> <chr> <chr> <dbl>
1 Afghanistan Men National 1985 20.2
2 Afghanistan Men Rural 1985 19.7
3 Afghanistan Men Urban 1985 22.4
4 Afghanistan Men National 2017 22.8
5 Afghanistan Men Rural 2017 22.5
6 Afghanistan Men Urban 2017 23.6Question opportunity: Why exactly are there 2,400 rows now?
Great! Our data is very usable now in this format!
Now that our data is clean and in a format that we can work with, we can start to take a look at the data and how different groups might compare.
Going back to one of our original questions,
Is there a difference between rural and urban BMI estimates around the world?
We want to use a statistical test to identify if there are differences in the average BMI estimates between e.g. rural and urban groups. First, let’s explore our data a bit.
National_BMI_1985 Rural_BMI_1985 Urban_BMI_1985 National_BMI_2017
Min. :18.20 Min. :17.7 Min. :19.30 Min. :20.20
1st Qu.:21.50 1st Qu.:20.9 1st Qu.:22.43 1st Qu.:24.07
Median :24.00 Median :23.4 Median :24.30 Median :26.15
Mean :23.68 Mean :23.3 Mean :24.24 Mean :26.00
3rd Qu.:25.23 3rd Qu.:25.1 3rd Qu.:25.40 3rd Qu.:27.43
Max. :34.30 Max. :34.1 Max. :34.40 Max. :35.30
NA's :6 NA's :2
Rural_BMI_2017 Urban_BMI_2017 Country Sex
Min. :19.80 Min. :21.50 Length:400 Length:400
1st Qu.:23.30 1st Qu.:24.80 Class :character Class :character
Median :26.00 Median :26.30 Mode :character Mode :character
Mean :25.61 Mean :26.37
3rd Qu.:27.30 3rd Qu.:27.60
Max. :35.00 Max. :35.40
NA's :8 NA's :2 We can see that the mean BMI estimates are larger in 2017. However, is this a meaningful increase? This is difficult to determine without a statistical test.
Lets look at the data stratified by Sex:
National_BMI_1985 Rural_BMI_1985 Urban_BMI_1985 National_BMI_2017
Min. :18.20 Min. :17.70 Min. :19.70 Min. :21.00
1st Qu.:22.00 1st Qu.:21.20 1st Qu.:22.85 1st Qu.:24.40
Median :24.25 Median :23.90 Median :24.60 Median :26.10
Mean :24.00 Mean :23.59 Mean :24.63 Mean :26.38
3rd Qu.:25.50 3rd Qu.:25.40 3rd Qu.:25.60 3rd Qu.:28.00
Max. :34.30 Max. :34.10 Max. :34.40 Max. :35.30
NA's :3 NA's :1
Rural_BMI_2017 Urban_BMI_2017 Country Sex
Min. :20.40 Min. :21.60 Length:200 Length:200
1st Qu.:23.70 1st Qu.:25.30 Class :character Class :character
Median :26.20 Median :26.40 Mode :character Mode :character
Mean :25.98 Mean :26.84
3rd Qu.:27.60 3rd Qu.:28.25
Max. :35.00 Max. :35.40
NA's :4 NA's :1 National_BMI_1985 Rural_BMI_1985 Urban_BMI_1985 National_BMI_2017
Min. :18.50 Min. :18.40 Min. :19.30 Min. :20.20
1st Qu.:21.10 1st Qu.:20.80 1st Qu.:22.10 1st Qu.:23.30
Median :23.60 Median :23.20 Median :24.10 Median :26.20
Mean :23.36 Mean :23.01 Mean :23.84 Mean :25.61
3rd Qu.:25.00 3rd Qu.:24.90 3rd Qu.:25.10 3rd Qu.:27.10
Max. :33.70 Max. :32.60 Max. :34.00 Max. :34.30
NA's :3 NA's :1
Rural_BMI_2017 Urban_BMI_2017 Country Sex
Min. :19.80 Min. :21.5 Length:200 Length:200
1st Qu.:22.60 1st Qu.:23.8 Class :character Class :character
Median :25.75 Median :26.3 Mode :character Mode :character
Mean :25.24 Mean :25.9
3rd Qu.:27.00 3rd Qu.:27.2
Max. :34.60 Max. :34.2
NA's :4 NA's :1 It looks like mean BMIs have increased in all regions for both men and women, but let’s look deeper with some statistical tests.
In order to apply a statistical test to compare the means, one of the first things to do is to explore the frequency of the different observed values. One way to summarize the frequency of different observed values is the frequency distribution, which can be shown in a table or a plot. See here for more information about distributions.
There are other types of distributions too, such as probability distributions, which describe the probability the occurrence of different possible outcomes. Probability distributions can be defined for discrete data (e.g. flipping a coin as heads or tails) and continuous data (e.g. BMI estimates). If we are talking about discrete data, the relative frequency of the different categories (heads or tails) can be used to define a discrete probability distribution. We can simply assign a probability to each category, e.g. \(P(heads) = 0.5\) and \(P(tails) = 0.5\).
Continuous probability distributions work in a similar way, except instead of assigning a probability to each category or type of observation (e.g. heads or tails), we assign a probability to an interval of continuous values (e.g. BMI between 22 and 26). This is because it’s not useful to construct a distribution that assigns a probability to each possible outcome (e.g. a BMI of 22.11, 22.12, etc) with such high precision. In this case, we would have to assign a very small probability to every possible BMI outcome.
The normal (or Gaussian) distribution is a type of theoretical continuous probability distribution that is frequently used to approximate the distribution of many naturally occurring measurements, such as that of BMI. We say that a random quantity \(X\) (e.g. BMI) is normally distributed with mean \(m\) and standard deviation \(s\) if the proportion of values in a given interval \((a,b)\) (e.g. BMI between 22 and 26) can be computed using this mathematical formula:
\[P(a < x < b) = \int_a^b \frac{1}{\sqrt{2\pi} s} e^{-\frac{1}{2}\big(\frac{x-m}{s}\big)^2} \]
Translating this to our case study, if our BMI data follow a normal distribution with mean of \(m = 24\) and standard deviation of \(s = 1\), then the distribution of different BMI outcomes should be centered around that mean value of 24 and would look something like this:
norm_BMI_ex_data <-
tibble(norm_data = rnorm(n = 200, mean = 24, sd = 1))
hist(pull(norm_BMI_ex_data, norm_data))
Here we are simulating 200 observations from a normal distribution with a mean of \(m = 24\) and standard deviation of \(s = 1\) using the rnorm() function.
Alternatively to using base R, we can plot the frequency of simulated observations using the geom_hist() function of the ggplot2 R package. The ggplot2 package creates plots by using layers. Notice in the following code how there is a plus sign between the ggplot() function and the geom_histogram() function. With ggplot2 we select what data we would like to plot using the first function (ggplot()) and then we add on additional layers of complexity (these layers can even involve different data).
Here, we are using the column called norm_data within norm_BMI_ex_data. As we noted, with ggplot2 we must first specify the data to plot with the ggplot() function. Then, the geom_* function specifies what type of plot to create (e.g. geom_histogram() create a histogram).
We could also make this same plot by piping the norm_BMI_ex_data data into the ggplot() function:
We will see later how we can add many layers to plots with ggplot2. For additional information on using ggplot2, see this case study.
Now instead of using simulated data from a normal distribution, let’s make a histogram of our BMI data:
Interesting! The data looks like it has one mode, but it does not look completely symmetric. In fact it looks like it might be what is called right skewed. See here for more information about skewed distributions.
To check we can calculate the mean and median of the BMI data using the base summarize() function of the dplyr package. If the mean and median are very similar, then we will know that the data isn’t skewed. Note that it is necessary to filter out the NA values first.
# A tibble: 1 x 2
mean median
<dbl> <dbl>
1 24.9 24.9OK, so they are essentially the same, mean= 24.87 and median = 24.9, so this does not fit the hallmark of right-skewed data where the the mean is greater than the median.
However so far we have made a histogram of the BMI data in the long format which combines BMI estimates for different subgroups of people, so rural women, rural men, urban women, urban men, etc. are all combined. If we want to make group comparisons, we need to look at the distributions for the specific groups that we are evaluating. So, we need to pick a particular population of individuals and graph their BMI, for example just the national estimates for women and for men in 1985.
The easiest way to do this is to use some other functions of the ggplot2 package, particularly the facet_wrap() function, which will allow us to look at differences in the distribution of the BMI estimates stratified by Year, Sex, and Region. We can sequentially divide our plots by deeper levels using multiple variables and the + plus sign within facet_wrap().
Some of these plots look pretty similar to a normal distribution, like the Women 2017 Urban plot (although it might be right-skewed, again we could check the mean and median to determine that). However some of the other plots, like the Men 2017 Urban plot, look very different. These plots show what is called a bimodal distribution, meaning that there appear to be two peaks (or modes), or in other words two different values (and values near each of these) that are the most frequent.
Another useful plot to create when exploring distributions is something called a “quantile-quantile” plot (or Q-Q plot for short). When we talk about a quantile, we are talking about dividing up the distribution of the data into roughly equal portions where roughly the same number of observations fall into each portion. For example, if you divide your data into 100 quantiles, you can think about this as percentiles, but you could also divide your data into 10 quantiles and these would be called deciles.
Why Q-Q plots? This plot allows us to compare the quantiles of two distributions together: (1) quantiles of a known theoretical distribution (like the normal distribution) compared to (2) quantiles of the distribution of our data. If the quantiles from these two distributions line up in the plot, then that is a visual piece of evidence that our data follow that theoretical distribution (like the normal distribution).
How does this work? To do this we will plot the quantiles of our data on the y-axis and the quantiles of the theoretical normal distribution on the x-axis. If the quantiles line up then we can say that our data is fairly normal. See here for more information about Q-Q plots.
Using the stat_qq() function of the ggplot2 package, we can easily create a Q-Q plot for our data randomly sampled from a normal distribution (“sample”) and compare it to the quantiles from a normal distribution (“theoretical”). The default comparison distribution for these functions is the normal distribution, so we don’t need to specify it in our code.
The stat_qq_line() function is used to add a line (computes the slope an intercept) on the plot. If the points lie on this straight line, this is evidence that the data have a normal distribution, not that the data have a particular normal distribution. Not surprisingly, we see that the points mostly fall on the line, indicating that the quantiles are fairly similar between the observed and theoretical data.
There are some points that are a bit further from the line as we get to the extreme quantiles.
Note that this plot is comparing our norm_BMI_ex_data which was randomly sampled from a normal distribution with mean 24 and standard deviation 1 to a standard normal distribution with mean 0 and standard deviation 1, which is why the values on the x-axis are -3 to 3 instead of 21 to 27.
As expected we see that about half the points are below our mean of 24 on the y-axis.
We can specify a specific mean and standard deviation for our theoretical distribution in our Q-Q plot; here we specify a normal distribution with mean 24 and standard deviation 1 (so that the theoretical distribution looks more like our data):
norm_BMI_ex_data %>%
ggplot(aes(sample = norm_data)) +
stat_qq(distribution = qnorm, dparams = c(mean = 24, sd = 1)) +
stat_qq_line(distribution = qnorm, dparams = c(mean = 24, sd = 1))
Notice that the plot looks the same except for the values on the x-axis have changed, meaning our interpretation will be the same – it looks like the values are consistent with a normal distribution.
Now, let’s take a look at our actual BMI data:
BMI_long %>%
ggplot(aes(sample =BMI)) +
stat_qq() +
stat_qq_line() +
facet_wrap(~ Sex + Year + Region)
This is a little hard to see, some plots suggest those subsets of the data are similar to a normal distribution, but others seem to be skewed. Let’s take a closer look.
For the sake of simplicity, we are going to focus on the subset of data for the women. This is because women were identified in the paper to have fastest increases in BMI over time. If we want to compare two subsets of women (e.g. 1985 versus 2017 or rural versus urban), then we need to look at the distribution of each of these subsets.
BMI_long %>%
filter(Sex == "Women") %>%
ggplot(aes(sample = BMI)) +
stat_qq() +
stat_qq_line() +
facet_wrap(~ Year + Region)
We can see that at the extremes of our quantiles, for most of our data, our tails are not very similar to the theoretical distribution. We can can also see that our data is right- or positive-skewed, because the extreme values are above the line.
The rural data looks more normal than the urban data. You might be asking, why do we care if the data look more normal or not? It turns out if you want to ask questions like “Is there a difference in the average BMI estimates from women in 1985 compared to women in 2017?”, we will need to use some of type of statistical test. The most commonly used statistical tests that can answer that question rely on an assumption about the data you have, namely that your data have a similar shape as a theoretical normal distribution. So, if we want to use those powerful statistical tests that rely on this assumption of “normality”, then we need to check that assumption before using those tests.
We can conclude from these plots that the assumption of normality appears to be violated in our data.
In the previous section, we hinted at the existence of statistical tests to be able to assess if there is a difference in the average BMI estimates between two groups (e.g. national BMI estimates in rural areas and national BMI estimates in urban areas). One such test is the Student’s \(t\)-test, which can be used to determine if two group means are different. However, as mentioned before, this test depends on certain assumption of our data.
In this section, we are going to talk about this test and other tests that do not require the same assumptions about our data.
Let’s remind ourselves of one of our original questions,
Is there a difference between rural and urban BMI estimates around the world?
In hypothesis testing, we are interested in comparing two different hypotheses: a “null” hypothesis (can be thought of like a baseline e.g. the means between two groups are the same) compared to an “alternative” hypothesis (e.g. the means between two groups are different). We are going to ask if there is enough evidence in our data to reject the null hypothesis.
Let’s try to formalize this a bit.
Suppose we want to test whether there is a difference between rural and urban BMI estimates using the data from 1985. In our case, then, we want to test whether the mean BMI estimate among the rural areas is the same as the mean BMI estimate in the urban areas. If we call the true unknown means of the two groups \(\mu_U\) and \(\mu_R\), for the urban and rural areas, respectively, then we can define the null hypothesis that there is no difference in the two means:
\[ H_0: \mu_U = \mu_R \]
In contrast, we also define an alternative hypothesis that there is a difference between the mean BMI estimates:
\[ H_a: \mu_U \neq \mu_R \]
The idea behind a hypothesis test is that we assume the null hypothesis is true and we use our data to help us identify if there is enough evidence to reject the null hypothesis.
This is similar to the idea of assuming that individuals are not guilty until proven otherwise. If there is not enough evidence in the data, then we say we “fail to reject the null hypothesis”.
You might be asking, “But we do not know what the true group means are?” That is correct. Instead, we must estimate these means based on the data that we have. For example, in our data, we can estimate the average BMI estimates for women in rural and urban regions in 1985. Here we will use the summarize() function from the dplyr package again, this time using the argument na.rm = TRUE to remove NA values.
BMI_long %>%
filter(Year == 1985, Sex == "Women", Region != "National") %>%
group_by(Region) %>%
summarize(avg_bmi = mean(BMI, na.rm = TRUE)) # A tibble: 2 x 2
Region avg_bmi
<chr> <dbl>
1 Rural 23.6
2 Urban 24.6Here we see that that in 1985, the average BMI estimate for women across urban areas is 24.6 compared to 23.6 across rural areas. It looks like BMI estimates are, on average, higher in urban areas compared to rural areas. However, now we want to use a statistical test to assess if this difference is statistically meaningful or a difference of this size could be found due to random chance.
The two-sample \(t\)-test is a common way to determine if the means of two groups are different. It does this by calculating a observed test-statistic, \(t\), that looks at the difference in the sample means of the two groups divided by an estimate of the variability in the difference in sample means. Here \(s_{pooled}\) is standard deviation of the data when pooling the samples together. \[ t = \frac{\bar{X}_U-\bar{X}_R}{s_{pooled}\sqrt{\frac{1}{n_U} + \frac{1}{n_R}}} \] If the null hypothesis is true, and the means for urban and rural areas are the same, then the sample means in the numerator will be similar and this test statistic will be close to 0. However, if the test statistic is not close to 0, this provides evidence the means for the urban and rural areas are difference, and so we may choose to reject the null hypothesis.
To formally decide whether to reject the null hypothesis, we compare the observed test statistic to the distribution of possible test statistics that could occur if the null hypothesis were true using what we call a p-value. The p-value tells us the probability of observing a test statistic as or more extreme than the one we observed if the null hypothesis is true. If the p-value is small, say less than 0.05, then our observed test statistic is pretty rare if the null hypothesis is true. This would lead us to reject the null hypothesis and conclude that the means of the two groups are different. We will talk more about p-values, and in particular our choice of 0.05 as a comparison for the p-value, a little later.
Here we perform a two-sample \(t\)-test comparing estimate for rural women in 1985 to estimates for urban women in 1985. We use the var.equal = TRUE option here to use the pooled standard deviation \(s_{pooled}\) from the formula above.
t.test(pull(filter(BMI_long, Sex == "Women",
Year == "1985",
Region == "Rural"), BMI),
pull(filter(BMI_long, Sex == "Women",
Year == "1985",
Region == "Urban"), BMI),
var.equal = TRUE)
Two Sample t-test
data: pull(filter(BMI_long, Sex == "Women", Year == "1985", Region == "Rural"), BMI) and pull(filter(BMI_long, Sex == "Women", Year == "1985", Region == "Urban"), BMI)
t = -3.8952, df = 394, p-value = 0.0001152
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-1.5744694 -0.5182378
sample estimates:
mean of x mean of y
23.58782 24.63417 In this case the p-value for the test is 0.00012, which is quite small. Here would would reject the null hypothesis that the means in the two groups are equal and conclude there are differences in mean BMI estimates for women in urban and rural areas.
But wait! We need to check the assumptions of the statistical test we used to be sure it is valid. The two-sample \(t\)-test relies on several assumptions:
If these assumptions are violated, this doesn’t necessarily mean we can’t perform a statistical test. It just means we may need to consider the following options:
R by using the var.equal = TRUE option in the t.test() function.Are these assumptions met in our case? We’ve already seen, using Q-Q plots, that these BMI estimates may not be normally distributed. We haven’t checked the assumption of equal variance yet. But it seems clear that the last assumption (of independent samples) is violated in our case! The BMI data for rural areas comes from the same countries and the BMI data for urban areas, meaning the samples are not independent and are instead what we call paired. We need to account for this because paired values may be more similar to one another because they come from the same country. Note that this would also be true for the male and female measurements from the same country or the values in the same countries from 1985 and later in 2017.
When we perform a paired \(t\)-test, instead of looking at the means of the two samples separately, we combine observations in a pair (from the same country in our case) to calculate a difference in values between the two samples. For example, we are really testing the bmi_diff column in the data below:
BMI_long %>%
filter(Year == 1985, Sex == "Women", Region != "National") %>%
group_by(Country) %>%
summarize(bmi_rural = BMI[Region == "Rural"],
bmi_urban = BMI[Region == "Urban"],
bmi_diff = bmi_rural - bmi_urban) # A tibble: 200 x 4
Country bmi_rural bmi_urban bmi_diff
<chr> <dbl> <dbl> <dbl>
1 Afghanistan 20.1 23.2 -3.10
2 Albania 26.1 25.9 0.2
3 Algeria 23.3 24.8 -1.5
4 American Samoa 34.1 34.4 -0.300
5 Andorra 25.4 25.2 0.200
6 Angola 20.9 22.7 -1.8
7 Antigua and Barbuda 24 25.4 -1.40
8 Argentina 23.3 24.5 -1.20
9 Armenia 25.2 25.5 -0.3
10 Australia 24.7 24.4 0.3
# … with 190 more rowsThis means our hypotheses are slightly different from the two-sample \(t\)-test. In this case we are testing the differences among the pairs of observations and how close these differences are to zero. (A mean difference of zero would mean the BMI estimates are the same in the rural and urban areas.) If we let \(\mu_d\) represent the true mean difference between the paired observations of the two groups, then our null hypothesis is that the mean of the differences is equal to zero: \[ H_0 : \mu_d = 0 \]
And the alternative hypothesis is that the mean of the differences is not equal to zero: \[ H_a : \mu_d \neq 0 \]
We can perform a paired \(t\)-test by adding the option paired = TRUE to the t.test() command. Note that we don’t need to specify anything about equality of variances now, because there is only one variance, the variance of the paired differences!
t.test(pull(filter(BMI_long, Sex == "Women",
Year == "1985",
Region == "Rural"), BMI),
pull(filter(BMI_long, Sex == "Women",
Year == "1985",
Region == "Urban"), BMI),
paired = TRUE)
Paired t-test
data: pull(filter(BMI_long, Sex == "Women", Year == "1985", Region == "Rural"), BMI) and pull(filter(BMI_long, Sex == "Women", Year == "1985", Region == "Urban"), BMI)
t = -14.095, df = 195, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-1.1870263 -0.8956268
sample estimates:
mean of the differences
-1.041327 Question opportunity: Looking at the t value, was global BMI lower in Rural or Urban areas in 1985?
Let’s come back to the other \(t\)-test assumptions now. With a paired \(t\)-test, we no longer need to assess the assumption about equal variance, since we are really testing the mean of one sample (of paired differences) rather than comparing two samples. However, we still have to check normality and independence of observations. We can rephrase these now to be about the paired differences rather than the two separate samples:
If we assume that measurements between different countries are independent, this second assumption is not violated. Therefore, we only now need to check normality of these paired differences.
We can start by making a Q-Q plot of the paired differences.
BMI_long %>%
filter(Year == 1985, Sex == "Women", Region != "National") %>%
group_by(Country) %>%
summarize(bmi_rural = BMI[Region == "Rural"],
bmi_urban = BMI[Region == "Urban"],
bmi_diff = bmi_rural - bmi_urban) %>%
ggplot(aes(sample = bmi_diff)) +
stat_qq() +
stat_qq_line()
We see that these paired differences do not appear to be normally distributed because the points do not lie on the straight line. We can also see this by making a histogram of these paired differences:
BMI_long %>%
filter(Year == 1985, Sex == "Women", Region != "National") %>%
group_by(Country) %>%
summarize(bmi_rural = BMI[Region == "Rural"],
bmi_urban = BMI[Region == "Urban"],
bmi_diff = bmi_rural - bmi_urban) %>%
ggplot(aes(x = bmi_diff)) +
geom_histogram()
This histogram appears to show show a bimodal distribution, meaning two groups of countries – most countries have a rural BMI estimates lower than urban BMI estimates (bmi_diff < 0), while a smaller group of countries have urban BMI estimates higher than urban BMI estimates (bmi_diff >0).
First we will try transform our data to make it more normally distributed. One way to do this is to take the logarithm of the data values. Logarithms make extreme values less extreme, so can be helpful when data has a right skew.
# create a new dataset with the paired differences
# create a log version of the difference variable using mutate()
BMI_diff_log <- BMI_long %>%
filter(Year == 1985, Sex == "Women", Region != "National") %>%
group_by(Country) %>%
summarize(bmi_rural = BMI[Region == "Rural"],
bmi_urban = BMI[Region == "Urban"],
bmi_diff = bmi_rural - bmi_urban) %>%
mutate(log_diff=log(bmi_diff))Notice the warning message that In log(bmi_diff) : NaNs produced. This is because some of our differences were less than zero, and the logarithm of a negative number is not defined. We can fix this issue by adding the same positive number, say 5, to each difference before taking the logarithm:
BMI_diff_log <- BMI_long %>%
filter(Year == 1985, Sex == "Women", Region != "National") %>%
group_by(Country) %>%
summarize(bmi_rural = BMI[Region == "Rural"],
bmi_urban = BMI[Region == "Urban"],
bmi_diff = bmi_rural - bmi_urban) %>%
mutate(log_diff=log(bmi_diff + 5))Let’s look at a Q-Q plot and histogram of this transformed data:
The transformed data still does not appear to be normally distributed and the bimodal nature of the distribution is still apparent in the histogram. What should we do? The \(t\)-test is fairly robust to violations of the normality assumption if the sample size is relatively large, due to what is called the central limit theorem, which states that as samples get larger, the sample mean has an approximate normal distribution.
With our sample size of 200, we could consider using the paired \(t\)-test even with our non-normal data. However, if we had a smaller sample, we might still need to consider another alternative since the log-transformed data still doesn’t look normal. We might want to consider a non-parametric test in that case.
When the normality assumption of a paired \(t\)-test is violated, we can consider a nonparametric alternative instead. What’s a non-parametric test? Parametric tests are based on assumptions about the distribution of the underlying data, while nonparametric tests do not rely on these distributional assumptions.
See here for more information about the difference between these two classes of tests.
The nonparametric alternative to the paired \(t\)-test is called the Wilcoxon signed-rank test. This test doesn’t require the same normality assumption as the \(t\)-test, so can be considered with the data does not appear to be normally distributed and particularly when the number of samples is low. The null hypothesis for this test is that the median of the paired differences is 0.
To do a Wilcoxon test, we use the wilcox.test() function:
wilcox.test(pull(filter(BMI_long, Sex == "Women",
Year == "1985",
Region == "Rural"), BMI),
pull(filter(BMI_long, Sex == "Women",
Year == "1985",
Region == "Urban"), BMI),
paired = TRUE)
Wilcoxon signed rank test with continuity correction
data: pull(filter(BMI_long, Sex == "Women", Year == "1985", Region == "Rural"), BMI) and pull(filter(BMI_long, Sex == "Women", Year == "1985", Region == "Urban"), BMI)
V = 1562, p-value < 2.2e-16
alternative hypothesis: true location shift is not equal to 0With this nonparametric test, we still see a small p-value, indicating that BMI estimates between urban and rural areas are different. We would conclude that there’s a difference in the median BMI estimates between these types of regions.
Let’s take a look back at the questions we wanted to answer: 1) Is there a difference between rural and urban BMI estimates around the world? In particular, how does this look for females? 2) How have BMI estimates changed from 1985 to 2017? Again, In particular, how does this look for females? 3) How do different countries compare for BMI estimates? In particular, how does the United States compare to the rest of the world?
It’s important to note that ultimately we wanted to test 4 different tests to answer question 1 and 2:
So far we have only looked at the first test. Let’s perform all 4 of these tests using the non-parametric Wilcoxon signed-rank test:
wilcox.test(pull(filter(BMI_long, Sex == "Women",
Year == "1985",
Region == "Rural"), BMI),
pull(filter(BMI_long, Sex == "Women",
Year == "1985",
Region == "Urban"), BMI),
paired = TRUE)
Wilcoxon signed rank test with continuity correction
data: pull(filter(BMI_long, Sex == "Women", Year == "1985", Region == "Rural"), BMI) and pull(filter(BMI_long, Sex == "Women", Year == "1985", Region == "Urban"), BMI)
V = 1562, p-value < 2.2e-16
alternative hypothesis: true location shift is not equal to 0wilcox.test(pull(filter(BMI_long, Sex == "Women",
Year == "2017",
Region == "Rural"), BMI),
pull(filter(BMI_long, Sex == "Women",
Year == "2017",
Region == "Urban"), BMI),
paired = TRUE)
Wilcoxon signed rank test with continuity correction
data: pull(filter(BMI_long, Sex == "Women", Year == "2017", Region == "Rural"), BMI) and pull(filter(BMI_long, Sex == "Women", Year == "2017", Region == "Urban"), BMI)
V = 2750.5, p-value < 2.2e-16
alternative hypothesis: true location shift is not equal to 0wilcox.test(pull(filter(BMI_long, Sex == "Women",
Year == "1985",
Region == "Rural"), BMI),
pull(filter(BMI_long, Sex == "Women",
Year == "2017",
Region == "Rural"), BMI),
paired = TRUE)
Wilcoxon signed rank test with continuity correction
data: pull(filter(BMI_long, Sex == "Women", Year == "1985", Region == "Rural"), BMI) and pull(filter(BMI_long, Sex == "Women", Year == "2017", Region == "Rural"), BMI)
V = 273, p-value < 2.2e-16
alternative hypothesis: true location shift is not equal to 0wilcox.test(pull(filter(BMI_long, Sex == "Women",
Year == "1985",
Region == "Urban"), BMI),
pull(filter(BMI_long, Sex == "Women",
Year == "2017",
Region == "Urban"), BMI),
paired = TRUE)
Wilcoxon signed rank test with continuity correction
data: pull(filter(BMI_long, Sex == "Women", Year == "1985", Region == "Urban"), BMI) and pull(filter(BMI_long, Sex == "Women", Year == "2017", Region == "Urban"), BMI)
V = 189, p-value < 2.2e-16
alternative hypothesis: true location shift is not equal to 0In each of these tests, the very small p-value that would lead us to conclude that there is a difference between the two groups we are comparing. However, when we test multiple questions like this, we need to correct for the multiple tests that we are performing.
The more we test data, the more likely we are to see a significant finding just by random chance even when the finding is in fact not real. It all comes down to probability.
You may recall that probability is the chance that an event will happen and it ranges from 0 to 1. You can also think of it as a percentage of 0% (probability =0) to 100% (probability = 1) of chance that an event will happen.
[source]
If we consider flipping a coin, the probability of an event happening (heads) is always equal to 1 ( or 100% chance) minus the probability of the event not happening (tails instead of heads) because the total probability can not exceed 100%.
probability of 1 (100% chance) = the probabilty of heads (50% chance)
+ the probability of tails (50% chance)
1 = P(heads) + P(tails)
probability of one outcome = 1 - the probability of the other outcome
P(heads) = 1 - P(tails)
or
P(tails) = 1 - P(heads)If we think about the probability of getting significant or nonsignificant results we can think of it like this:
1 = P(signifiant results)
+ P(no significant results)
P(significant results) = 1 - P(no significant results)Alpha is the amount of risk that you are willing to accept that a statistical test result is actually a false positive, in other words we reject the null hypothesis when the null hypothesis is true. This is also called type 1 error. A significance threshold of .05 is customary and it means that we are accepting a 5% chance that our result is actually a false positive.
alpha = probability of false significant results
or
α = P(Making a type 1 error)
Remember that in our statistical tests if the p-value < alpha it is considered significant and if p-value > alpha it is considered not significant
But what really is the p-value?
The p in p-value stands for probability - the probability that we would obtain the statistics (for example the t in our student \(t\)-tests based on the means of our groups of comparison) just by random chance alone. Therefore a p-value of 0.02 means that there is a 2 percent chance that out data may look the way it does just because of random chance and not because there is really a difference in the means of the groups of interest.
So then if alpha is the threshold for the p-value for obtaining false significant results then the probability of not making incorrect conclusions is 1 - alpha:
P(Not making an error) = 1 - α[1] 0.95OK so if we use an alpha of .05 we are accepting that 95% of the time we will not make a Type 1 error and 5 % of the time we will just by random chance.
Taking this one step further:
P(Making an error) = 1 - P(Not making an error)
P(Making an error) = 1 - (1 - α)
Here we can see that this checks out:
[1] 0.05So what happens if we perform multiple tests?
The probability of not making a type 1 error would remain the same for each test. Therefore we need to multiply the probabilities together each time to determine the over all probability of making an error across multiple tests. See here about why we multiply probabilities together.
P(Not making an error in m tests) = (1 - α)^m
P(Making at least 1 error in m tests) = 1 - (1 - α)^mLet’s consider if we performed 10 different statistical tests and if we as usual considered the significance threshold alpha of .05:
So the probability of getting 1 significant result with 10 tests is:
#P(Making at least 1 error in 10 tests) = 1 - ((1 - α)(1 - α)(1 - α)(1 - α)(1 - α)(1 - α)(1 - α)(1 - α)(1 - α)(1 - α))
#this is the same as:
#P(Making at least 1 error in 10 tests) = 1 - (1 - α)^10
1- (1-.05)^10[1] 0.4012631So there is a 40% chance that that there will be a significant finding simply due to random chance alone.
What about 100 tests?
[1] 0.9940795Yikes!! That is almost a 100% chance that there will be a significant finding simply due to chance alone!
Much of this explanation is described in this lecture.
One way to correct for this is multiple testing issue is using the Bonferroni method.
In this method we would divide our significance threshold (generally 0.05) by the number of tests.
[1] 0.0125Our new significance threshold is now 0.0125. Thus our p-values should be less than this value for us to reject the null that there is no difference in means. In our 4 Wilcoxon tests, our p-values were less than 0.0125. So we see a significant difference in the means of the groups after multiple testing correction for our four different tests.
Again we will utilize ggplot2 to create plots to look at the directional and magnitude of the differences in BMI we are interested in. If you need additional information please see here. The top two lines of the code for the following plots, filter the data to only specific values of interest. Then we layer what is called a jitter on top of a box plot. A jitter is essentially a dot plot but with some variation on the location of the points so that they do not line up vertically which can make the individual points difficult to see.
Our first main question was: Is there a difference between rural and urban BMI estimates around the world? Let’s look at the national mean BMI estimates for each of the years:
BMI_long %>%
filter(Sex == "Women",
Year == "1985",
Region != "National") %>%
ggplot(aes(x = Region, y = BMI)) +
geom_boxplot() +
geom_jitter(width = .3) 
# the width determines how wide the points are plotted
BMI_long %>%
filter(Sex == "Women",
Year == "2017",
Region != "National") %>%
ggplot(aes(x = Region, y = BMI)) +
geom_boxplot() +
geom_jitter(width = .3)
Our 2nd main question was: How have BMI estimates changed from 1985 to 2017? Let’s look at the change in rural and urban mean BMI estimates over time:
#this time we will add titles for Rural and Urban using ggtitle()
BMI_long %>%
filter(Sex == "Women",
Year %in% c("1985", "2017"),
Region == "Rural") %>%
# filtering this way allows us to keep data from both years
ggplot(aes(x = Year, y = BMI)) +
geom_boxplot() +
geom_jitter(width = .3) +
ggtitle("Rural")
BMI_long %>%
filter(Sex == "Women",
Year %in% c("1985", "2017"),
Region == "Urban") %>%
ggplot(aes(x = Year, y = BMI)) +
geom_boxplot() +
geom_jitter(width = .3) +
ggtitle("Urban")
Let’s put the plots together to see how the change over the years differs between the regions. We will use again use facet_wrap() to do this:
BMI_long %>%
filter(Sex == "Women",
Year %in% c("1985", "2017"),
Region %in% c("Rural", "Urban")) %>%
# filtering this way allows us to keep data from both years
ggplot(aes(x = Year, y = BMI)) +
geom_boxplot() +
geom_jitter(width = .3) +
facet_wrap(~ Region)
Indeed the change in BMI over time looks bigger in the rural areas!
One third main questions was: How do the different countries compare? Or in other words what do the individual dots represent in our box plots? We will take a look using geom_label():
BMI_long %>%
filter(Sex == "Women",
Year %in% c("1985", "2017"),
Region == "Rural") %>%
ggplot(aes(x = Year, y = BMI)) +
geom_boxplot() +
geom_jitter(width = .3) +
geom_label(aes(label = Country))
If we include all country names this is a bit too much… so perhaps we should focus on just the extreme BMI values using filter() function of thedplyr package. We will also use the ggrepel package to have our labels not overlap each other.
Part of the third question was: How does the United States of America compare? So let’s label the data points for the United states.
BMI_long %>%
filter(Sex == "Women",
Year %in% c("1985", "2017"),
Region == "Rural") %>%
ggplot(aes(x = Year, y = BMI)) +
geom_boxplot() +
geom_jitter(data=BMI_long %>%
filter(Sex == "Women",
Year %in% c("1985", "2017"),
Region == "Rural") %>%
filter(BMI>31 | BMI<19 | Country == "United States of America"),
aes(x =Year, y =BMI),
width = .02) +
ggrepel::geom_text_repel(data=BMI_long %>%
filter(Sex == "Women",
Year %in% c("1985", "2017"),
Region == "Rural") %>%
filter(BMI>31 | BMI<19 | Country == "United States of America"),
aes(x =Year, y =BMI, label = Country))
And let’s fill the box plots with color and outline in black:
BMI_long %>%
filter(Sex == "Women",
Year %in% c("1985", "2017"),
Region == "Rural") %>%
ggplot(aes(x = Year, y = BMI)) +
geom_boxplot(outlier.shape = NA, color = "black", aes(fill = Year)) +
geom_jitter(data=BMI_long %>%
filter(Sex == "Women",
Year %in% c("1985", "2017"),
Region == "Rural") %>%
subset(BMI>31 | BMI<19 | Country == "United States of America"),
aes(x =Year, y =BMI),
width = .02) +
ggrepel::geom_text_repel(data=BMI_long %>%
filter(Sex == "Women",
Year %in% c("1985", "2017"),
Region == "Rural") %>%
subset(BMI>31 | BMI<19 | Country == "United States of America"),
aes(x =Year, y =BMI, label = Country))
Let’s take a look at all the data together, including the data for men:
ggplot(BMI_long, aes(x = Year, y = BMI, col = Region)) +
geom_boxplot(outlier.shape = NA, color = "black" , aes(fill = Region)) +
facet_grid(~ Sex) 
That’s useful, but let’s look at the individual points and include our country labels, to do so lets change our United States of America label to USA:
BMI_long[["Country"]] <-BMI_long[["Country"]] %>%
str_replace( pattern = "United States of America", replacement = "USA")
ggplot(BMI_long, aes(x = Year, y = BMI, col = Year)) +
geom_boxplot(outlier.shape = NA, color = "black" , aes(fill = Sex)) +
facet_grid(~ Sex + Region) +
geom_hline(yintercept=30, linetype="dashed", color = "red", size =1) +
# This will add a horizontal dashed line to indicate the obesity BMI threshold
geom_jitter(data=subset(BMI_long),
aes(x =Year, y =BMI),
width = .2, size =2, shape =21, color = "black", fill = "gray") +
# This will add the individual country data points
# The shape 21 allows for a different fill and outline color
# The width determines how wide the jitter points are plotted
geom_jitter(data=subset(BMI_long, Country == "USA"),
aes(x =Year, y =BMI),
width = .02, size =12, shape =21, color = "black", fill = "gray") +
# This will add points that are larger for the USA data
geom_text(data=subset(BMI_long, Country == "USA"),
aes(x =Year, y =BMI,label=Country),
color = "black") +
# This will add USA labels to the USA points
theme(legend.position = "none",
# This is useful for removing the legend
axis.text.x = element_text(size = 15,angle = 30),
# this changes the size and angle of the x axis point labels
axis.text.y = element_text(size = 20),
axis.title.y = element_text(size =15),
axis.title.x = element_text(size =15),
strip.text.x = element_text(size = 15)) +
# This changes the size of x axis labels for the facet
ggtitle( "Differences in BMI Over Time and Across Region Type and Gender") +
# Add a plot title
scale_fill_manual(values=c("dodgerblue", "orchid2"))
Here we can see that overall BMI appears to be increasing globally over time. Additionally we can see that this is occurring not just in urban areas, but also in rural areas. The US is consistently above the median in all strata of the data. In general, the female data shows higher BMI values than the male data. The rural USA BMI appears to be higher than the urban BMI for both men and women. Many countries have average BMI estimates above the obesity threshold of 30. Thus it appears that education and outreach programs for weight management should focus on both rural and urban areas and both genders. Education and assistance for women may be especially helpful.
How does the rate of change in BMI differ between groups? Which group might especially need attention?
First let’s calculate the differences in BMI from 2017 and 1985 and add this to our BMI_long data object. This will create new variables called Rural_difference, Urban_difference, and National_difference.
BMI_numeric[["Rural_difference"]] <-
BMI_numeric[["Rural_BMI_2017"]] - BMI_numeric[["Rural_BMI_1985"]]
BMI_numeric[["Urban_difference"]] <-
BMI_numeric[["Urban_BMI_2017"]] - BMI_numeric[["Urban_BMI_1985"]]
BMI_numeric[["National_difference"]] <-
BMI_numeric[["National_BMI_2017"]] - BMI_numeric[["National_BMI_1985"]]
BMI_diff_long <- BMI_numeric %>%
select(Country:National_difference) %>%
gather(key = Type, value = Difference , Rural_difference:National_difference)
BMI_diff_long <- BMI_numeric %>%
select(Country: National_difference) %>%
pivot_longer( names_to = "Difference" , cols = Rural_difference:National_difference)
head(BMI_diff_long)# A tibble: 6 x 4
Country Sex Difference value
<chr> <chr> <chr> <dbl>
1 Afghanistan Men Rural_difference 2.8
2 Afghanistan Men Urban_difference 1.2
3 Afghanistan Men National_difference 2.6
4 Afghanistan Women Rural_difference 3.5
5 Afghanistan Women Urban_difference 3.1
6 Afghanistan Women National_difference 3.80Let’s replace “United states of America” with “USA” and make a plot with this data to compare the change in BMI:
BMI_diff_long[["Country"]] <-BMI_diff_long[["Country"]] %>%
str_replace( pattern = "United States of America", replacement = "USA")
ggplot(BMI_diff_long, aes(x = Difference, y = value, col = Difference)) +
geom_boxplot(outlier.shape = NA, color = "black" , aes(fill = Difference)) +
facet_grid(~ Sex) +
geom_jitter(data=subset(BMI_diff_long),
aes(x =Difference, y =value),
width = .2, size = 2, shape = 21, color = "black", fill = "gray") +
# This will add the individual country data points
# The shape 21 allows for a different fill and outline color
# The width determines how wide the jitter points are plotted
geom_jitter(data=subset(BMI_diff_long, Country == "USA"),
aes(x = Difference, y = value),
width = .02, size = 12, shape = 21, color = "black", fill = "gray") +
# This will add points that are larger for the USA data
geom_text(data=subset(BMI_diff_long, Country == "USA"),
aes(x =Difference, y = value, label = Country), color = "black") +
# This will add USA labels to the USA points
theme(legend.position = "none",
# This is useful for removing the legend
axis.text.x = element_text(size = 15,angle = 30),
# this changes the size and angle of the x axis point labels
axis.text.y = element_text(size = 20),
axis.title.y = element_text(size = 15),
axis.title.x = element_text(size = 15),
strip.text.x = element_text(size = 15)) +
# this changes the size of x axis labels for the facet
ggtitle( "Change in BMI Over Time and Across Region Type and Gender")
We can now see that the rate of change from 1985 to 2017 appears to be larger in the women compared to the men in all regions. The group with the largest increase in the USA is the women living in rural areas.
Let’s check the difference with some statistical tests. Remember that these rural differences are measured on the same countries as the urban differences, so we will need to use paired tests to make these comparisons. We can either consider the paired \(t\)-test or the Wilcoxon signed-rank test. To decide between them, we’ll look at the normality of the paired differences in change over time.
Let’s compare the change over time in rural communities to the change in urban communities among women and separately among men.
# lets look at the normality of these paired differences
# first we need to calculate the paired differences of change over time
BMI_diff_long %>%
filter(Difference != "National_difference") %>%
group_by(Country, Sex) %>%
summarize(change_rural = value[Difference == "Rural_difference"],
change_urban = value[Difference == "Urban_difference"],
change_diff = change_rural - change_urban) %>%
ggplot(aes(x=change_diff)) +
geom_histogram() +
facet_wrap(~ Sex)
# interesting...the spread of the difference in change over time is much
# broader for women than for men
# for women, the differences look reasonably normal, spread between -1 and +1
# for men there appear to be two groups of countries; one with difference in change over time
# centered around 0 and the other with difference in change over time between 1 and 2.
# Let's look at Q-Q plots:
BMI_diff_long %>%
filter(Difference != "National_difference") %>%
group_by(Country, Sex) %>%
summarize(change_rural = value[Difference == "Rural_difference"],
change_urban = value[Difference == "Urban_difference"],
change_diff = change_rural - change_urban) %>%
ggplot(aes(sample = change_diff)) +
stat_qq() + # to add the line we need to also include stat_qq_line()
stat_qq_line() +
facet_wrap(~ Sex)
# the distribution of differences in change over time does not appear normal
# for either men or womenBased on the non-normality of these paired differences, we will use a Wilcoxon signed-rank test to compare changes over time between urban and rural areas separately for men and women:
wilcox.test(pull(filter(BMI_diff_long,
Sex == "Women",
Difference == "Rural_difference"), value),
pull(filter(BMI_diff_long,
Sex == "Women",
Difference == "Urban_difference"), value),
paired = TRUE)
Wilcoxon signed rank test with continuity correction
data: pull(filter(BMI_diff_long, Sex == "Women", Difference == "Rural_difference"), value) and pull(filter(BMI_diff_long, Sex == "Women", Difference == "Urban_difference"), value)
V = 11304, p-value = 0.00327
alternative hypothesis: true location shift is not equal to 0wilcox.test(pull(filter(BMI_diff_long,
Sex == "Men",
Difference == "Rural_difference"), value),
pull(filter(BMI_diff_long,
Sex == "Men",
Difference == "Urban_difference"), value),
paired = TRUE)
Wilcoxon signed rank test with continuity correction
data: pull(filter(BMI_diff_long, Sex == "Men", Difference == "Rural_difference"), value) and pull(filter(BMI_diff_long, Sex == "Men", Difference == "Urban_difference"), value)
V = 9159, p-value = 0.08256
alternative hypothesis: true location shift is not equal to 0Using a p-value threshold of 0.05, we would conclude there’s a difference in change over time between urban and rural areas for women (p = 0.00327) but not for men (p = 0.08256).
However, now we have performed six comparisons (4 earlier) so we should apply our multiple testing correction and use a new p-value threshold of:
[1] 0.008333333With the threshold adjusted for multiple comparisons, we would still conclude there’s a difference in the change over time between urban and rural areas for women but not for men.
Interestingly, although we don’t see a difference in change over time for men, our histogram of the differences between urban and rural areas for men highlighted two groups of countries – one centered around zero (no difference) and one above 1 (showing a difference):
BMI_diff_long %>%
filter(Difference != "National_difference") %>%
group_by(Country, Sex) %>%
summarize(change_rural = value[Difference == "Rural_difference"],
change_urban = value[Difference == "Urban_difference"],
change_diff = change_rural - change_urban) %>%
ggplot(aes(x=change_diff)) +
geom_histogram() +
facet_wrap(~ Sex)
In fact for these countries, the difference between rural and urban areas for men is higher than the difference for women for almost all of the countries!
We can look to see which countries have such a large difference in change over time between rural and urban areas for men:
BMI_diff_long %>%
filter(Sex == "Men", Difference != "National_difference") %>%
group_by(Country) %>%
summarize(change_rural = value[Difference == "Rural_difference"],
change_urban = value[Difference == "Urban_difference"],
change_diff = change_rural - change_urban) %>%
filter(change_diff > 1)# A tibble: 21 x 4
Country change_rural change_urban change_diff
<chr> <dbl> <dbl> <dbl>
1 Afghanistan 2.8 1.2 1.60
2 American Samoa 2 0.2 1.80
3 Bangladesh 2.3 0.600 1.7
4 Bhutan 2.8 1.10 1.7
5 Cook Islands 3.10 1.4 1.70
6 Fiji 3 1.60 1.4
7 French Polynesia 2.5 0.8 1.70
8 India 2.6 0.7 1.90
9 Kiribati 2.10 0.5 1.60
10 Marshall Islands 2.5 0.900 1.6
# … with 11 more rowsFor women, we did see a difference in the change over time between urban and rural areas. There appear to be specific countries where BMI shows a particular increase especially for women. Which countries are those? How does that compare with the US? Clearly the US is among the countries with the highest differences.
# A tibble: 6 x 4
Country Sex Difference value
<chr> <chr> <chr> <dbl>
1 USA Men Rural_difference 3.10
2 USA Men Urban_difference 3.1
3 USA Men National_difference 3.10
4 USA Women Rural_difference 3.5
5 USA Women Urban_difference 3.40
6 USA Women National_difference 3.4 In the US, focus should be placed on both urban and rural women to improve this public health issue.
Here we can see the countries that have the largest differences in BMI from 1985-2017:
# A tibble: 6 x 4
Country Sex Difference value
<chr> <chr> <chr> <dbl>
1 Egypt Women Rural_difference 5.9
2 Honduras Women National_difference 5.60
3 Honduras Women Rural_difference 5.5
4 Egypt Women National_difference 5.10
5 Saint Lucia Women Rural_difference 5
6 Honduras Women Urban_difference 4.9 However it is important to see what the mean BMI values were for these countries in 2017. It could be that the average was underweight in 1985… let’s take a look.
BMI_long %>%
filter(Country %in% pull(filter(BMI_diff_long, value > 4.9), Country),
Sex %in% pull(filter(BMI_diff_long, value > 4.9), Sex),
Year =="2017") %>%
arrange(- BMI)# A tibble: 9 x 5
Country Sex Region Year BMI
<chr> <chr> <chr> <chr> <dbl>
1 Egypt Women Urban 2017 32.3
2 Egypt Women National 2017 31.7
3 Egypt Women Rural 2017 31.3
4 Saint Lucia Women National 2017 30.5
5 Saint Lucia Women Rural 2017 30.5
6 Saint Lucia Women Urban 2017 30.5
7 Honduras Women Urban 2017 28.3
8 Honduras Women National 2017 27.7
9 Honduras Women Rural 2017 26.9#Here is the data just for the changes in BMI in Egypt
filter(BMI_diff_long, Country == "Egypt") %>%
arrange(- value)# A tibble: 6 x 4
Country Sex Difference value
<chr> <chr> <chr> <dbl>
1 Egypt Women Rural_difference 5.9
2 Egypt Women National_difference 5.10
3 Egypt Women Urban_difference 4.20
4 Egypt Men Urban_difference 2.60
5 Egypt Men Rural_difference 2.5
6 Egypt Men National_difference 2.5 Thus rural women in Egypt had the greatest increase in BMI from 1985 to 2017 in this data (5.9) and the average BMI is now over the obesity threshold of 30.
The data suggests that rural women in Egypt and other countries may especially benefit from dietary resources and our interventions and programs to assist with weight management.
Now let’s make a plot that summarizes our findings. To do this we will simplify the other plots we made and then combine them together.
# simplified national means plot
Means_plot<-BMI_long %>%
filter(Sex %in% c("Men", "Women"),
Year %in% c("1985", "2017"),
Region == "National") %>%
ggplot(aes(x = Sex, y = BMI)) +
geom_boxplot(outlier.shape = NA, color = "black" , aes(fill = Sex)) +
scale_fill_manual(values=c("dodgerblue", "orchid2")) +
facet_grid(~ Year) +
geom_hline(yintercept=30, linetype="dashed", color = "red", size =1) +
geom_jitter(data=BMI_long %>%
filter(Sex %in% c("Men", "Women"),
Year %in% c("1985", "2017"),
Region == "National"),
aes(x =Sex, y =BMI),
width = .2, size =2, shape =21,
color = "black", fill = "gray") +
geom_jitter(data=subset(BMI_long %>%
filter(Sex %in% c("Men", "Women"),
Year %in% c("1985", "2017"),
Region == "National",
Country == "USA")),
aes(x =Sex, y =BMI),
width = .02, size =12, shape =21,
color = "black", fill = "gray") +
# This will add points that are larger for the USA data
geom_text(data=subset(BMI_long%>%
filter(Sex %in% c("Men", "Women"),
Year %in% c("1985", "2017"),
Region == "National", Country == "USA")),
aes(x =Sex, y =BMI,label=Country),
color = "black") +
# This will add USA labels to the USA points
theme_linedraw() +
# This will make thebackground of the plot white and the facet labels black
theme(legend.position = "none",
# This is useful for removing the legend
axis.text.x = element_text(size = 15,angle = 30, vjust = 0.5),
# this changes the size and angle of the x axis point labels
axis.text.y = element_text(size = 20),
axis.title.y = element_text(size =15),
axis.title.x = element_text(size =15),
strip.text.x = element_text(size = 15))
#Simplified difference plot
BMI_diff_long[["Difference"]] <-BMI_diff_long[["Difference"]]%>%
str_replace( pattern = "_difference", replacement = "")
Diff_plot<-BMI_diff_long %>%
filter(Sex %in% c("Men", "Women"),
Difference != "National") %>%
ggplot( aes(x = Difference, y = value, col = Difference)) +
geom_boxplot(outlier.shape = NA, color = "black" , aes(fill = Sex)) +
scale_fill_manual(values=c("dodgerblue", "orchid2")) +
facet_grid(~ Sex) +
geom_jitter(data=BMI_diff_long %>%
filter(Sex %in% c("Men", "Women"),
Difference != "National"),
aes(x =Difference, y = value),
width = .2, size = 2, shape = 21,
color = "black", fill = "gray") +
geom_jitter(data=subset(BMI_diff_long %>%
filter(Sex %in% c("Men", "Women"),
Difference != "National",
Country == "USA")),
aes(x = Difference, y = value),
width = .02, size =12, shape =21,
color = "black", fill = "gray") +
# This will add points that are larger for the USA data
geom_text(data=subset(BMI_diff_long %>%
filter(Sex %in% c("Men", "Women"),
Difference != "National"),
Country == "USA"),
aes(x = Difference, y = value,label=Country),
color = "black") +
# This will add USA labels to the USA points
theme_linedraw() +
# This will make the background of the plot white and the facet labels black
theme(legend.position = "none",
# This is useful for removing the legend
axis.text.x = element_text(size = 15,angle = 30, vjust = 0.5),
# this changes the size and angle of the x axis point labels
axis.text.y = element_text(size = 20),
axis.title.y = element_text(size =15),
axis.title.x = element_text(size =15),
strip.text.x = element_text(size = 15))
# this changes the size of x axis labels for the facet
#add labels
Diff_plot<-Diff_plot +
labs(title = "Change in BMI by region",
x = "",
y = "Change in BMI \n (1985 to 2017)") +
theme(title = element_text (size = 12, face = "bold"))
Means_plot <-Means_plot +
labs(title = "Mean National BMI over time",
x = "",
y = "Mean BMI") +
theme(title = element_text (size = 12, face = "bold"))
obesity_text<-tibble(Year=c(1985),BMI=c(31),Sex=c("Men"),label=c("Obesity"))
Means_plot <-Means_plot +
geom_text(data = obesity_text,
label=pull(obesity_text,label),
color = "red",
aes( fontface ="bold.italic", size = 13))
cowplot::plot_grid(Means_plot, Diff_plot, labels = c("A", "B"))
We could make a similar plot with the patchwork package:
Means_plot +labs(tag = 'A') + Diff_plot+ labs(tag = 'B') +
plot_layout(guides = 'collect', ncol = 2)
Great, now we have put two plots together using the plot_grid() function of the cowplot package. This way we can clearly communicate two messages. The first being that BMI has increased over time globally and that many countries including the United States of America are approaching a mean BMI that is above the obesity threshold of 30. We can also see that women on average have larger BMI values than males, but that both genders show increased levels over time. In the second plot we can see that the increase in BMI is not just happening in urban communities, but in both rural and urban communities and particularly in women.
Our plot visually explores all of our main questions: 1) Is there a difference between rural and urban BMI estimates around the world? In particular, how does this look for females? 2) How have BMI estimates changed from 1985 to 2017? Again, In particular, how does this look for females? 3) How do different countries compare for BMI estimates? In particular, how does the United States compare to the rest of the world?
Now let’s save our plot as a png file using the png() function:
png(here::here("img", "main_plot.png"), height = 1500, width = 2300, res = 300)
cowplot::plot_grid(Means_plot, Diff_plot, labels = c("A", "B"))
dev.off()quartz_off_screen
2 You might think that these two plots look a bit odd together, that if you look at the change in the boxplots in part A, that the data for Men shows a bigger change, however, we can’t tell from this plot how the individual countries are changing. Part B of the figure looks at the specific change for each country.
We could check our plot using the following code:
BMI_diff_long %>%
group_by(Sex, Difference) %>%
filter(!is.na(value))%>%
filter(Country == "USA") %>%
summarize(mean = mean(value), median = median(value))# A tibble: 6 x 4
# Groups: Sex [2]
Sex Difference mean median
<chr> <chr> <dbl> <dbl>
1 Men National 3.10 3.10
2 Men Rural 3.10 3.10
3 Men Urban 3.1 3.1
4 Women National 3.4 3.4
5 Women Rural 3.5 3.5
6 Women Urban 3.40 3.40BMI_diff_long %>%
group_by(Sex, Difference) %>%
filter(!is.na(value))%>%
summarize(mean = mean(value), median = median(value))# A tibble: 6 x 4
# Groups: Sex [2]
Sex Difference mean median
<chr> <chr> <dbl> <dbl>
1 Men National 2.26 2.2
2 Men Rural 2.23 2.30
3 Men Urban 2.05 2.10
4 Women National 2.37 2.70
5 Women Rural 2.38 2.5
6 Women Urban 2.21 2.60BMI_long %>%
group_by(Sex, Year, Region) %>%
filter(!is.na(BMI))%>%
filter(Country == "USA") %>%
summarize(mean = mean(BMI), median = median(BMI))# A tibble: 12 x 5
# Groups: Sex, Year [4]
Sex Year Region mean median
<chr> <chr> <chr> <dbl> <dbl>
1 Men 1985 National 25.8 25.8
2 Men 1985 Rural 26.1 26.1
3 Men 1985 Urban 25.7 25.7
4 Men 2017 National 28.9 28.9
5 Men 2017 Rural 29.2 29.2
6 Men 2017 Urban 28.8 28.8
7 Women 1985 National 25.7 25.7
8 Women 1985 Rural 26 26
9 Women 1985 Urban 25.6 25.6
10 Women 2017 National 29.1 29.1
11 Women 2017 Rural 29.5 29.5
12 Women 2017 Urban 29 29 BMI_long %>%
group_by(Sex, Year, Region) %>%
filter(!is.na(BMI))%>%
summarize(mean = mean(BMI), median = median(BMI))# A tibble: 12 x 5
# Groups: Sex, Year [4]
Sex Year Region mean median
<chr> <chr> <chr> <dbl> <dbl>
1 Men 1985 National 23.4 23.6
2 Men 1985 Rural 23.0 23.2
3 Men 1985 Urban 23.8 24.1
4 Men 2017 National 25.6 26.2
5 Men 2017 Rural 25.2 25.8
6 Men 2017 Urban 25.9 26.3
7 Women 1985 National 24.0 24.2
8 Women 1985 Rural 23.6 23.9
9 Women 1985 Urban 24.6 24.6
10 Women 2017 National 26.4 26.1
11 Women 2017 Rural 26.0 26.2
12 Women 2017 Urban 26.8 26.4BMI_numeric %>%
group_by(Sex)%>%
summarize(Nat_mean1985 = mean(National_BMI_1985), Nat_mean2017 = mean(National_BMI_2017))# A tibble: 2 x 3
Sex Nat_mean1985 Nat_mean2017
<chr> <dbl> <dbl>
1 Men 23.4 25.6
2 Women 24.0 26.4BMI_numeric %>%
group_by(Sex)%>%
summarize(Nat_med1985 = median(National_BMI_1985), Nat_med2017 = median(National_BMI_2017))# A tibble: 2 x 3
Sex Nat_med1985 Nat_med2017
<chr> <dbl> <dbl>
1 Men 23.6 26.2
2 Women 24.2 26.1We can see that the largest changes occur for Women.
Scroll through the output!
BMI_diff_long %>%
filter(value>3.5) %>%
arrange(-value) %>%
# this allows us to show the full output
print(n = 1e3)# A tibble: 117 x 4
Country Sex Difference value
<chr> <chr> <chr> <dbl>
1 Egypt Women Rural 5.9
2 Honduras Women National 5.60
3 Honduras Women Rural 5.5
4 Egypt Women National 5.10
5 Saint Lucia Women Rural 5
6 Honduras Women Urban 4.9
7 El Salvador Women Rural 4.8
8 Mexico Women Rural 4.8
9 Uzbekistan Women Rural 4.8
10 Dominican Republic Women Rural 4.80
11 Comoros Women Urban 4.70
12 Nicaragua Women Rural 4.70
13 Saint Lucia Women National 4.70
14 El Salvador Women National 4.6
15 Haiti Women National 4.6
16 Kyrgyzstan Women Rural 4.60
17 Costa Rica Women Rural 4.5
18 Guatemala Women Rural 4.5
19 Jamaica Women Rural 4.5
20 Qatar Women Rural 4.4
21 Bahrain Women Rural 4.40
22 Dominican Republic Women National 4.40
23 Guatemala Women National 4.40
24 Malaysia Women Rural 4.40
25 Nicaragua Women National 4.40
26 Uzbekistan Women National 4.40
27 Bangladesh Women National 4.3
28 Gabon Women National 4.3
29 Iraq Women Rural 4.3
30 Mexico Women National 4.3
31 Saint Kitts and Nevis Women Rural 4.3
32 Belize Women Rural 4.20
33 Comoros Women National 4.20
34 Costa Rica Women National 4.20
35 Kyrgyzstan Women National 4.20
36 Panama Women Rural 4.20
37 Tajikistan Women Rural 4.20
38 Tanzania Women Urban 4.20
39 Egypt Women Urban 4.20
40 Brunei Darussalam Women Rural 4.1
41 Gambia Women National 4.1
42 Ghana Women National 4.1
43 Haiti Women Rural 4.1
44 Nepal Women National 4.1
45 Bolivia Women National 4.10
46 Chile Women Rural 4.10
47 Grenada Women Rural 4.10
48 Saint Lucia Men Rural 4.10
49 Samoa Women Rural 4.
50 Antigua and Barbuda Women Rural 4
51 Bahamas Women Rural 4
52 Bangladesh Women Rural 4
53 Bolivia Women Rural 4
54 Comoros Women Rural 4
55 Dominica Women Rural 4
56 Jamaica Women National 4
57 Malaysia Women National 4
58 Mexico Women Urban 4
59 Nicaragua Women Urban 4
60 Panama Women National 4
61 Puerto Rico Women Rural 4
62 Saint Kitts and Nevis Women National 4
63 Saint Lucia Men National 4
64 Saint Lucia Women Urban 4
65 Suriname Women Rural 4
66 Turkmenistan Women Rural 4
67 Bhutan Women National 3.9
68 Gabon Women Urban 3.9
69 Ghana Women Urban 3.9
70 Mali Women Urban 3.9
71 Samoa Women Urban 3.9
72 El Salvador Women Urban 3.90
73 Georgia Women Rural 3.90
74 Guyana Women Rural 3.90
75 Libya Women Rural 3.90
76 Saint Vincent and the Grenadines Women Rural 3.90
77 Samoa Women National 3.90
78 Grenada Women National 3.8
79 Nepal Women Rural 3.8
80 Pakistan Women National 3.8
81 Yemen Women National 3.8
82 Afghanistan Women National 3.80
83 Brunei Darussalam Women National 3.80
84 Guatemala Women Urban 3.80
85 Paraguay Women Rural 3.80
86 Tajikistan Women National 3.80
87 Dominica Women National 3.7
88 Guyana Women National 3.7
89 Algeria Women Rural 3.70
90 Argentina Women Rural 3.70
91 Belize Women National 3.70
92 Benin Women Urban 3.70
93 Benin Women National 3.70
94 Bolivia Women Urban 3.70
95 Dominican Republic Women Urban 3.70
96 Gambia Women Urban 3.70
97 Iraq Women National 3.70
98 Mongolia Women Rural 3.70
99 Occupied Palestinian Territory Women Rural 3.70
100 Pakistan Women Rural 3.70
101 Saint Lucia Men Urban 3.70
102 Saint Vincent and the Grenadines Women National 3.70
103 Syrian Arab Republic Women Rural 3.70
104 Uzbekistan Women Urban 3.70
105 Yemen Women Rural 3.70
106 Antigua and Barbuda Women National 3.6
107 Costa Rica Women Urban 3.6
108 Liberia Women Urban 3.6
109 Malaysia Women Urban 3.6
110 Mali Women National 3.6
111 Papua New Guinea Women Rural 3.6
112 Togo Women Urban 3.6
113 Turkmenistan Women National 3.6
114 Uruguay Women Rural 3.6
115 Ecuador Women Rural 3.60
116 Malawi Women Urban 3.60
117 Qatar Women National 3.60We have evaluated BMI average estimates from 200 different countries around the world. To do so we imported data from a pdf using the pdftools package. We used tidyverse packages such as dplyr, stringr, and tidy to clean the data and get it in a workable format. Our statistical analysis focused on evaluating differences in BMI in females around the world across time and between rural and urban areas. We found a significant difference both between years and between the type of community among women globally using paired \(t\)-tests and nonparametric tests. Thus BMI has increased in women since 1985. Although BMI estimates are significantly higher in Urban areas compared to rural areas, BMI estimates have increased in both regions. We learned that there are assumptions and considerations to keep in mind when performing tests that compare two groups.
We learned that the two sample \(t\)-test relies on:
We learned that we can evaluate if our data is normally distributed by plotting the distribution and by creating Q-Q plots.
We learned that if our data is not normally distributed, we can consider these options:
We learned that if our data has unequal variance, we can consider using the Welch’s \(t\)-test.
Importantly, we learned that when we have paired data (are observations are not independent) and we intend to compare means, we should perform a paired test (paired \(t\)-test or the nonparametric Wilcoxon signed rank test. Our data used in this case study was paired because observations were taken for the same countries for different categories of populations - thus, we wanted to compare these populations within the same country (male vs female, rural vs urban etc.). Other examples would be in a case-control study (if samples are matched for various demographics) or a study with repeated measures (for example, symptom measures from the same individual before and after a treatment). We learned that paired \(t\)-test tests a different hypothesis the two-sample \(t\)-test. In this case we are testing the differences among the pairs of observations and how close these differences are to zero, thus there is only one variance, the variance of the paired differences and therefore, we do not need to worry about the equality of variance assumption like we do with the two-sample \(t\)test. We learned that we can perform a paired \(t\)-test by adding the option paired = TRUE to the t.test() command.
The paired \(t\)-test assumes that:
We learned that if the normality assumption is violated, we can use the Wilcoxon signed rank test, or if our data is large enough, due to the central limit theorem, the paired \(t\)-test should still be robust to deviations from normality.
Using the ggplot2 package we were able to visualize trends in the data. Importantly, the largest increase appears to be in the rural areas particularly for women. We were also able to identify how the US compared to other countries and which countries showed the largest increase in BMI over time. Analyses like this are important for defining which groups could benefit the most from interventions, education, and policy changes when attempting to mitigate public health challenges. You can see in the article that this data came from that many additional considerations would be involved to perform such an analysis.
Students can evaluate the change in BMI over time using the global data available for each year between 2015 and 2017. The data can be found here.
BMI
Long and Wide Data Formats
Distributions Normal Distribution Skewed Distributions Bimodal Distribution ggplot2
Q-Q Plots
Student \(t\)-test
Paired Data
Welch’s \(t\)-test
Parametric and Nonparametric Methods
Variance
Balanced Study Design
Independent Observations
Transformation
Permutation/Resampling Methods
Central Limit Theorem
Mood’s Two-Sample Scale Test
Wilcoxon Signed Rank Test
Wilcoxon Rank Sum Test
Two-sample Kolmogorov-Smirnov Test
Type 1 Error
p-value
Multiple Testing
Bonferroni Method of Multiple Testing Correction